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1. 


A  MATHEMATICAL  MODEL  OF  WATER  ENTRY.  (U/U) 


PRECIS 


1.  u  A  computer  simulation  of  the  v&ter  entry  of  an  axisymmetric  tody  with  or 
without  a  cruciform  tail  and  with  or  without  a  parachute  delivery  system  is 
described.  The  predictions  of  the  simulation  are  shown  to  agree  with  experi¬ 
mental  observations  of  water  entry  motion.  The  Fortran  program  which  implements 
this  model  is  listed. 


A 


CONCLUSIONS 


2.  Over  the  range  of  impact  velocities  (20  to  10  m/sec)  which  were  experi¬ 
mentally  investigated  using  a  full  scale  dummy  torpedo,  the  simulation  gave 
reasonable  agreement  with  the  measured  water  entry  behaviour.  In  addition  a 
limited  comparison  at  higher  impact  velocities  indicated  that  the  possibility 
of  applying  the  simulation  to  a  wider  range  of  entry  velocities  was  promising. 


3.  The  principal  areas  in  which  future  work,  could  be  carried  out  are  in 
improving  the  model  of  the  splash,  the  cavity  and  the  interaction  of  both  the 
body  and  the  parachute  with  the  cavity  flow  field  particularly  at  low  Froude 
numbers. 
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3. 


INTRODUCTION 


U .  During  the  passage  of  a  body  from  air  into  water  a  complex  series  of 
forces  act  upon  the  body.  The  experimentally  observed  phenomena  associated  with 
this  water  entry  are  described  in  detail  by  Waugh  and  Stubstad  (Ref  l).  These 
various  phases  of  water  entry  are  summarised  in  figure  1. 

5.  At  initial  impact  and  during  the  subsequent  flow  formation  momentum  is 
rapidly  transferred  from  the  body  to  the  water  in  order  to  create  a  velocity 
field  in  the  water.  During  initial  impact  very  high  local  pressures  are 
experienced,  the  amplitudes  of  which  are  dependent  upon  the  physical  properties, 
including  compressibility,  of  the  body,  of  the  water,  and  of  the  air.  However 
the  actual  force  and  moment  impulses  are  comparatively  independent  of  com¬ 
pressibility  effects. 

6.  The  flow  field  contains  regions  of  low  pressure  and  consequently  cavita¬ 
tion  occurs.  Initially,  figure  lb,  the  cavity  is  open  to  the  atmosphere. 

During  this  open  cavity  phase  the  pressure  in  the  cavity  under  the  body  may  be 
significantly  lower  than  the,  near  atmospheric,  pressure  in  the  upper  part  of 
the  cavity.  As  the  cavitation  number  increases  the  cavity  closes,  figure  lc , 
and  progressively  decreases  in  size. 

7.  Some  of  the  air  which  was  originally  sucked  into  the  cavity  is  lost  by 
entrainment  at  the  rear  of  the  cavity.  Depending  upon  the  initial  entry  con¬ 
ditions  of  the  body,  the  tail  may  momentarily,  intermittently,  or  continuously 
contact  the  wall  of  the  cavity,  figure  Id.  Finally  the  cavity  collapses,  the 
body  slips  out  of  any  remaining  air  bubble,  figure  le,  and  becomes  fully  wet. 

8.  The  work  which  is  described  here  seeks  to  predict  the  forces  which  are 
applied  to  t'ne  body  during  these  various  phases  and  thus  to  predict  the  result¬ 
ing  motion  of  a  body  during  water  entry.  Although  the  mathematical  model  is 

of  an  analytical  nature  many  of  the  coefficients  and  functional  relationships 
have  been  determined  empirically  by  comparison  with  the  work  of  other  experi¬ 
menters  and  by  comparison  with  a  series  of  measurements  which  were  made  speci¬ 
fically  to  assist  in  forming  and  validating  this  model. 

9.  It  is  intended  that  this  mathematical  model  should  be  applicable  to  the 
water  entry  of  any  axisymmetric  body,  with  or  without  an  axisymmetric  parachute. 
However  the  actual  simulation  was  derived  with  particular  reference  to  a  light¬ 
weight  torpedo.  A  torpedo  may  be  launched  from  a  surface  ship  using  above 
water  torpedo  tubes  and  a  typical  torpedo  is  shown  in  figure  2.  The  torpedo  may 
also  be  delivered  by  an  aircraft  using  a  parachute  as  shown  in  figure  3.  In 
addition  a  torpedo  may  be  fitted  with  a  frangible  nose  cap  which  is  designed 

to  attenuate  the  initial  shock  at  water  impact. 

THE  MATHEMATICAL  MODEL 


Reference  frames 


10.  The  orientation  of  a  body  in  space  is  usually  defined  by  the  three  Euler 
angles,  roll  ($),  pitch  (0),  and  yaw  (ifi),  however  this,  representation  possesses 
a  singularity  when  the  body  is  pitched  at  ninety  degrees.  The  quaternion  four 
parameter  system,  first  described  by  the  Irish  mathematician  Sir  W  R  Hamilton 
(Ref  2),  which  also  may  be  used  to  define  the  attitude  of  a  body  in  space 
overcomes  this  singularity. 

11,  If  a  moving  so  or  rignt,  handed  axes  (x,y,z),  fixed  in  a  bouy,  ale 

obtained  from  the  m  co-ordinate  sixes  (x  ,y  ,x„),  fixed  in  space,  by 

S  S  o 
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rotating  the  space  frame  of  reference  through  the  angle  y  about  the  unit 
vector  (a  ,  (3  )  then  the  quaternion  parameters,  .e,  describing  the  oriei 

tion  of  the  body  axes  relative  to  the  space  axes  are  defined  as:- 


e  =  cos  -r 

o  2 


e,  »  a  sm  ^ 
1  s  2 


e2  =  Ps  3in  2 


e  =  V  sin  — 

3  s  2 

It  may  be  seen  that:- 

P  -  1  (2) 

12.  In  terms  of  the  three  Euler  angles  the  quaternion  parameters  are  given 

t,-;- 

<t>  0  i>  ,  ■  6  .  6  .  \p 

e  =•  cos  -±  cos  -r  cos  +  Sin  -Jr  sm  -r  sin  3C 
o  2  2  2  2  2  2 

i  0  b  d)  .  0  .  \p 

e±  ~  sin  cos  ~  cos  ^  ~  cos  ^  sin  sin  ~ 

(3) 

<j>  .  6  ip  .  <b  0  .  \b 

e„  --  cos  -i  sm  —  cos  -r  +  sm  -v  cos  -e-  sm  -r 

d  d  d  d  d  d  d 

<i  0  .  ip  .  d>  .  f)  d) 

C3  "  GOR  2  003  2  ■':in  2  “  sin  2  sin  2  cos  O 

13.  The  transformation  matrix,  which  relates  the  body  frame  of  reference  to 
i.-i/i  space  frame  of  reference,  is: 


2  2  2 
e  +e  -e  -e 
0  12  3 


2cm  +2e  e 
3  0  12 


^1_-;  oe2 


2e, e  -2e  e„ 

12  o  3 

2  2  2  2 
e  -e  +e  -e 
o  I  2  3 


2eoV2e2e3 


2e  e_+2e,  e„ 
o  2  1  o 


2e,,e,-2e  e. 
2  3  o  1 


2  2  2  2 
e  -e,  -e0  +e 
0  12  3 


ic  terminology :  ■ 


T  .  "I 

L~’  ■*’ 


1  1  1 

X2  Y2  Z2 

x3  Y3  Z3_. 
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5. 


will  "be  used  for  convenience.  In  the  space  frame  of  reference  the  unit  vector 
X  is  the  tody  x  axis,  the  unit  vector  Y_  is  the  "body  y  axis  and  the  unit  vector 
Z  is  the  "body  z  axis. 

The  equations  of  motion 

15 .  The  origin  of  the  body  fixed  axes  is  chosen  so  that  the  x  axis  is  the  axis 
of  symmetry  of  the  body  and  so  that  the  co-ordinates  of  the  position  of  the 
centre  of  gravity  are  (0syg,2g).  If  the  body,  of  mass  m,  axial  moment  of 
inertia  lx  and  transverse  moment  of  inertia  Iy,  is  moving  with  linear  velocity 
(u,v,w)  and  angular  velocity  (p,q,r)  under  the  influence  of  external  forces 
(Fx,Fy,Fz)  and  moments  (Lx,Ly,Lz)  in  a  gravitational  field  g  then  the  equations 
of  motion  are:- 


m  ^u  -  vr  +  wq  +  yg(pq  -  r)  +  zg(pr  +  4)  -  gO 

m(v-vp  +  ur-  yg(r  +  p  )  +  zg(qr  -  p)  -  gY3j 

(2  2  \ 

w  -  uq  +  vp  +  yg(rq  +  p)  -  zg(p  +  q  )  -  gZ^ 

IXP  +  m  (yg(w  -  uq  +  vp)  -  2g(v  -  vp  +  ur)  +  g(zgY3  -  ygZ3)^ 
I„q  +  (l„  -  rp  +  m  za  (u  -  vr  +  wq  -  gX^) 


_y"  '“x  _y '  _jr  8 

-  )  pq  +  m  y_ 
y  x  ^ 

The  geometry  of  water  entry 


Izr  +  (Iy  -  Ix)  pq  +  m  yg  (vr  -  u  -  wq  +  gX3) 


*  F 


=  F 


=  F 


=  L 


=  L 


=  L 


(6) 


16.  In  order  to  predict  the  external  forces  and  moments  applied  to  the  body 
it  is  first  necessary  to  determine  which  parts  of  the  body  are  in  contact  with 
water . 


17.  The  shape  of  the  axially  symmetric  body  is  defined  by  a  table  of  axial 
distances,  x,  and  the  corresponding  radii,  R,  so  that  the  body  is  divided  into 
a  number  of  segments,  one  of  which  is  shown  in  figure  From  this  table  of 
values  the  following  segment  parameters  are  de fined 


segment 

'x'  co-ordinate 

=  x  = 

(xn 

segment 

radius 

=  R  = 

(Rn 

segment 

width 

=  dx  = 

x  - 
n 

'  Xn+1 

segment 

angle 

=  a  = 

tan 

«\.l  - 

(7) 


-  Rn)/dx) 


18.  The  geometrical  problem  which  must  be  solved  in  order  to  determine  which 
areas  of  the  body  surface  are  wet  is  shown  in  figure  5*  For  an  area  of  the 
body  to  be  in  contact  with  water  it  is  necessary  that  the  area  be  both  beneath 
the  sea  surface  and  not  in  a  region  of  cavitation. 


Sea  surface  condition 


19-  If  the  unit  vector  in  space  out  of  the  sea  surface  is  n^,  then  in  the  body 
co-ordinate  system  this  is;- 
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The  body  segment,  in  body  co-ordinates  relative  to  the  centre  of  the  segment 
may  be  defined  as  (o,  R  cos  g,  R  sin  g).  Hence  the  intersection  of  the 
segment  with  the  sea  surface  is  defined  by  the  two  roots  g^,  g^  of:- 

(0,  R  cos  g,  R  sin  g)  .  n  =  d  (9) 

where  d  is  the  normal  distance  between  the  plane  of  the  sea  surface  and  the 
centre  of  the  segment.  Equation  9  has  a  solution  if:- 


d 


<  1 


(10) 


When  there  is  no  solution  the  sign  of  d  determines  whether  the  segment  is  com¬ 
pletely  above  or  completely  below  the  sea  surface. 

Cavitation  condition 

2C .  The  pressure  distribution,  cavity  detachment  angles,  and  cavity  shapes, 
described  by  Knapp,  Daily  and  Hammitt  (Ref  3),  for  a  number  of  axisymmetric 
bodies  were  studied.  The  empirical  conditions  for  cavity  detachment  are  chosen 
to  be  that  the  cavitation  number,  a,  is  less  than  0.3  and  that  the  local  angle 
of  incidence,  i,  of  the  body  surface  to  the  flow  is:- 


iin(i)  -  0.152  sin  (i  )  +  O.608  sin^  (i  )  -  0.34  exp  (-  R  cos  a  ~) 

HI  IU  {IX 


dcu 


-  0.30a  +  0.28  sin2  (i  ) 

D 


(11) 


who r«.  i  is  the  maximum  local  angle  of  incidence  and  i,  is  the  incidence  of  the 
11:  b 

wn ol'.  body. 


21.  If  the  velocity  is  assumed  to  be  constant  over  each  segment  of  the  body 
then  the  local  incidence,  i^ ,  of  a  point  at  angular  position  g  on  a  segment  is:- 


,  2.  2.  2\2 


sir.  (i.,)  ~  (u  .iina  +  (v,cosg  +  w  sing)  cos  a)/(uc  tv^+w  _) 

I*'*  S  S  S  S  S  S 


(12) 


"vr (u  ,v  )  is  the  velocity  of  the  centre  of  the  segment.  By  equating  i6 
to  i  from  equation  (il)  the  cuvitating  sector  is  found  in  terms  of  g.  However 


3 1  :• 


u  sin  a  -  (v‘"+w  )  cos  a  >  sin  (i)  (u'+v  +w  ) 
s  s  s  s  s 


(13) 


then  there  is  no  cavitation  and  if:- 


i  1 

p  P  2  pops 

u  sin  a  +  (v'+w  )  cos  a  <  sin  (i)  (u  +v‘”+w‘") 
s  s  s  s  s  s 


(14) 


then  the  whole  segment  is  cavitating  and  a  cavity  is  assumed  to  be  thrown  off 
from  the  body  at  this  segment. 

22.  This  cavity  is  assumed  to  be  an  ellipsoid  of  revolution  with  its  major 
axis  aligned  to  the  flow  direction  and  with : - 
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7. 


semi  major  axis  =  a  =  0.4  sin(i)  (l  +  R  cos  a  ~)/o 


semi  minor  axis  =  b  =  R  ^  +  0.13a 
where  the  effective  segment  radius  R  f  is 

2  l 

R ~  R(u2+v2+w2)  /  (u  -  (v2+w2)  tan  a) 

S  S  S  S  S  S 


(15) 

(16) 


(17) 


The  origin  of  this  ellipsoid,  is  found  by  fitting  the  radius  of  the  cavity  to  the 
effective  radius  of  the  segment.  Initially,  after  impact,  the  size  of  this 
cavity  and  the  pressure  of  the  air  entrained  within  it  are  defined  as  empirical 
functions  of  the  distance  travelled  in  water. 

23.  If  an  ellipsoidal  cavity  is  present  then  the  wetted  sectors  of  segments 
which  are  downstream  of  the  inception  of  this  cavity  are  determined  by  solving 
for  the  intersection  of  each  body  segment  with  the  cavity. 

The  external  forces  and  moments 

24.  The  external  forces  on  the  body  are  divided  into  those  forces  resulting 
from  quasi  steady  state  pressures  and  those  forces  which  are  associated  with 
'added  mass'  effects. 


Steady  state  uress^ire  forces 

25.  The  mass  of  water,  of  density  p,  displaced  by  each  wetted  sector  is:- 

dm  =  pR2dx  (02  -  -  sin  ({Jg-f »1))/2  (l8) 

therefore,  in  the  body  frame  of  reference,  the  buoyancy  force  on  each  segment 
is :  — 


dF  =  dm  g  (X3,  Y3,  Z3) 


(19) 


26.  The  dynamic  pressure  on  the  surface  of  the  body  is  defined  to  be  propor¬ 
tional  to  the  square  of  the  normal  component  of  the  velocity  of  the  surface. 

The  normal  component  of  velocity,  V0 ,  at  angular  position  8  on  a  segment  is  the 

p 

scalar  product  of  the  velocity  of  the  surface  and  the  unit  vector  normal  to  the 
surface  at  that  position  ie:- 

VQ  =  ((u,v,w)  +  ((p,q,r)  x  (x,Rcos0 ,Rsin$) ) )  .  (sina,  cosRcosa,  sinftcosa) 

(20) 

The  magnitude  of  the  pressure  force,  dF,  on  tne  elemental  area,  RdB  dx/cos  a, 
is 


dF 


Vg  Cn 

p  D  cos  a 


vrV>  r»y*  r* 


n 


the  force  coefficient. 


(21) 
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27*  If  a  constant  velocity,  V  ,  is  chosen  to  be  of  the  same  order  as  V  . 

6o  3 


then  to  a  first  approximation:- 


v;  ■  \  <2vs  -  v  (22) 

2 

Substituting  this  linear  expression  for  V  into  equation  (2l),  resolving  dF 

p 

into  the  body  frame  of  reference  and  integrating  over  the  whole  segment  yields :• 

e2 

dpx  =  I  p  CD  R  ^  tan  °  f  <2Vg  ■  v6  >  d3 


C. 

p  CD  Vg  R  dx  /  (2Vg  -  Vg  )  cos  3 


d?z  =  t  P  °D  V£?  E  ^  /  <2Vg  “  Vg  >  sin  E  a* 

o  j  0 

31 

The  moment,  dL ,  of  dF  about  the  origin  is:- 

oL  =  (x,  R  cos  6,  R  sin  6)  x  dF  (2M 

it:  dL  ~  ( h  tan  a  -  x)  (0,  dF  ,  -dF  )  (2?) 

—  z  y 

The  external  forces  and  moments  applied  to  the  body  by  quasi  steady  state 
pressures  are  obtained  by  summing  dF  and  dL  over  all  of  the  segments  of  the 

V:-dy . 

A  j d eci  muss  fore e s 

jo.  The-  importance  of  forces  associated  with  the  momentum  of  the  flow  field 
surrounding  the  partially  wet  body  during  water  entry  was  first  emphasised  by 
von  Karnian  (Ref  U).  The  linear  momentum  M,  and  the  angular  momentum,  H,  of  the 
flow  field  surrounding  the  body  may  be  defined  by:- 


M 

I:  n. 

__  z_J  [_  u 


X. 

X , 

X. 

X. 

X. 

X. 

u 

V 

w 

p 

q 

r 

Y. 

Y, 

Y, 

Y. 

Y, 

Y. 

u 

V 

w 

P 

q 

r 

Z. 

Z. 

Z, 

Z. 

z. 

Z. 

u 

V 

V 

P 

q 

r 

K. 

K. 

K. 

K. 

K. 

K. 

u 

V 

V 

P 

q 

r 

M. 

M. 

M. 

M. 

M. 

M. 

*vr 

V 

P 

q 

i 

N. 

N. 

N, 

N. 

N. 

N. 

_  u 

V 

V 

P 

q 

r_ 
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9. 


The  six  by  six  array  is  the  added  mass  matrix  which  is  discussed  in  detail  by 
Imlay  (Ref  5).  The  resulting  forces  and  moments  acting  on  a  moving  body  are:- 


dM 


1  “  +  ((p»q»r)  X  M) 


(27) 


+  ((u,v,w)  x  M)  +  ( (p,q,r)  x  H)  (28) 

to  evaluate  the  added  mass  matrix  it  is  assumed  that  associated 
with  every  element  of  area,  Rdgdx/cos  a,  is  a  volume  of  fluid,  Rhdgdx/cos  a, 
so  that  h  is  a  measure  of  the  distribution  of  added  mass  upon  the  body.  In  addi¬ 
tion  it  is  assumed  that  only  the  component  of  velocity  normal  to  the  surface 
imparts  momentum  to  the  associated  volume  of  water.  The  element  of  linear 
momentum  normal  to  the  surface  at  angular  position  8  on  a  segment  is  therefore 


dH 


-  =  dt 


and 

29 •  In  order 


dM  =  V„  p  Rh  d8  dx/cos  a 
p 


(29) 


where  V  is  given  by  equation  (20).  Integrating  over  a  whole  segment  gives  the 
p 

fluid  linear  momentum  associated  with  the  vetted  sector  of  the  segment  expressed 
in  the  body  frame  of  reference 

g2 

dM  =  p  R  h  dx  tan  a  f  V  d$ 
x  J  8 

e. 


dM  =  p  E  h  dx  T  V.  cos  8  d8 
y  J  p 

ex 

B2 

dM  =  p  R  h  dx  f  V.  sin  8  d8 
z  J  p 


(30) 


'.ngular  momentum  about  the  body  origin  is:- 
dH  =  (x,  R  cos  8,  R  sin  6)  x  dM  (31) 

ie  dH  =  (R  tan  a  -  x)  (0,  dM^,  -dM.^.)  (32) 

30.  The  coefficients  of  u,v,w,p,q  and  r  in  equations  (30)  and  (32)  are 
equated  to  the  identical  coefficients,  comprising  the  added  mass  matrix,  in 
equation  (26)  and  hence  the  contributions  of  each  wetted  sector  to  all  36  of  the 
added  mass  derivatives  are  obtained.  The  total  added  mass  matrix  is  formed  by 
summing  the  individual  contributions  from  each  segment  of  the  body. 


The  solution  of  the  equations  of  motion 
dM  dH 

31.  Noting  that  —  and  •—  contain  the  rates  of  change  of  the  added  masses  in 
dt  at 

addition  to  acceleration  terms,  the  total  forces  and  moments  both  due  to  the 
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added,  mass  effects,  equations  (27)  and  (28),  and  due  to  the  steady  state 
pressures,  equations  (19),  (23)  and  (25),  are  summed  linearly  and  substituted 
into  the  equations  of  motion,  equation  (6),  which  may  then  be  expressed  in  the 
form : - 


[•]  ::  ■[■]*•[>] 


where  A  is  a  six  by  six  matrix  defining  the  total  inertias  of  the  system,  B 
1 a  six  by  one  matrix  defining  the  steady  state  external  forces  and  C  is  a 
s.ix  by  one  matrix,  defining  the  changes  of  momentum  which  have  occurred  during 
the  rime  interval,  dt.  The  equations  of  motion  are  expressed  with  dt  as  a 
multiplicand  in  order  that  the  stable  numerical  integration  of  these  equations 
may  be  performed  through  the  indeterminate  accelerations  associated  with  the 
wetting  of  an  incompressible  body  by  an  incompressible  liquid. 

.12 .  At  each  cycle  of  the  numerical  integration  of  the  equations  of  motion 
the  six  simul  taneous  linear  equations  (33)  are  solved  to  yio;ld  the  velocity 
increments  from  which  the  linear  and  angular  velocities  of  the  body  are  up¬ 
date  i.  In  addition  the  orientation  and  position  in  space  (x  ,yo,z  )  of  the 
boo $  arc  obtained  by  integrating  the  kinematic  relationships 


-q  -r“] 


r  -q 


i  y.  I  =  T a  I  S  V  S 


Pi  ?c  as  si  on  of  the  y."  '.meric  a],  simulation 

33.  A  t'0 STRAP  program  which  generates  and  integrates  equations  3.3,  3^  and 
33  is  listed  in  the  appendix, 

3I1 ,  This  program  requires  the  shape  of  the  body  and  the  added  mass  distribu¬ 
tion  to  be  supplied  as  input  data.  The  added  mass  distribution  is  somewhat 
subjective,  however  if  ip  helpful  to  consider  some  examples  v/hich  will  assist 
in  estimating  the  added  mass  distribution  parameter,  h.  Lamb  (Ref  6,  p  Ihk) 
indicates  that  on  each  side  of  a  flat  disc  the  added  mass  distribution  is 
given  by:- 
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i 


O  Q  o  2 

h  *  |  (R  -  R||)  (36) 

where  R  is  the  radius  of  the  disc  and  R^  is  the  radial  position  at  which  h  is 
defined.  On  p.155  of  the  same  reference  it  is  shown  that : - 

h  *  R/2  (37) 

for  a  sphere  of  radius  R  and  that 

h  =  R  (38) 

for  the  sides  of  a  long  cylinder  of  radius  R. 

35.  In  addition  to  the  basic  model  described  in  the  previous  paragraphs  the 
simulation  Also  includes  the  effects  of  a  horizontal  steady  wind  and  a  simple 
sea  motion,  it  is  assumed  that  the  body  is  small  compared  to  velocity  gradients 
in  the  sea  and  that  the  sea  velocity  potential,  <J>,  may  be  represented  by 

t  ^  n  s  \ 

U\ei 

$  *  ac  exp  (-wzq)  cos  w  (xo~ct)  (39) 

where  a  is  the  wave  amplitude,  c  is  the  wave  celerity,  and  w  is  the  wave  fre¬ 
quency.  Typical  values  of  a,  c  and  v  are  tabulated  by  Lofft  and  Price  (Ref  7). 

36.  The  model  allows  the  addition  of  cruciform  tail  fins  and/or  a  shroud 

ring  tail.  These  tail  surfaces  ire  assumed  to  have  a  linear  relationship  between 
force  and  incidence  at  small  angles  of  incidence. 

37.  The  underpressure  effect  which  occurs  during  the  initial  period  of  oblique 
water  entry  is  represented  in  the  simulation  by  a  local  pressure  reduction  in 
the  cavity  on  the  underside  of  the  nose  of  the  body. 

38.  It  was  experimentally  observed  that  the  tail,  when  in  contact  with  the 
cavity  wall,  experiences  an  upward  force  which  is  thought  to  be  due  to  the 
gravity  effect  described  by  Knapp,  Daily  and  Hammitt  (Ref  3,  p.25l).  This 
effect  is  represented  in  the  simulation  by  an  additional  upward  velocity  field 
superimposed  on  the  rear  of  the  cavity  at  low  Foude  numbers. 

39.  When  the  body  is  fully  surrounded  by  a  single  fluid,  before  impact  or 
after  deep  cavity  collapse,  then  the  external  forces  on  the  body  are  represented 
by  force  derivatives  in  the  usual  way  (eg  Ref  8,  p.196). 

40 .  The  simulation  has  a  provision  for  the  shape  of  the  body  to  change  after 
a  prescribed  amount  of  kinetic  energy  has  been  dissipated  during  water  entry. 

This  facility  may  be  used  to  represent  a  frangible  nose  cap. 

1*1.  The  simulation,  as  listed,  will  not  be  valid  for  body  angles  of  incidence 
greater  than  about  135°  as  any  axial  cavity  from  the  tail  will  not  be  modelled 
correctly. 

1(2.  The  mathematical  model  described  on  the  previous  pages  may  be  applied  to 
an  ax i symmetric  parachute.  The  simulation  listed  in  the  appendix  allows  a 
parachute  to  be  attached  to  the  body  via  elastic  rigging  lines. 

1(3.  The  additional  rorces  applied  by  the  rigging  lines  to  both  the  parachute 
and  the  body  are  obtained  by  deriving  the  strain  and  strain  rate  of  each  line 
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from  the  known  position  and  velocity  of  the  body  relative  to  the  parachute. 
The  equations  of  motion  of  the  body  and  of  the  parachute  are  then  evaluated 
independently. 


14 .  In  the  previous  paragraphs  the  principles  of  the  water  entry  simulation 
were  described.  The  reader  who  wishes  to  explore  the  detailed  implementation 
of  these  principles  may  do  so  by  studying  the  appendix. 


THE  EXPERIMENTAL  MEASUREMENTS 


15.  The  difficulties  associated  with  scaling  water  entry  behaviour  are  dis¬ 
cussed  by  Knapp,  Daily  and  Hammitt  (Ref  3,  p.548).  In  view  of  the  many 
uncertainties  associated  with  extrapolating  small  scale  model  measurements 

up  to  full  scale  the  experimental  measurements  required  to  improve  and  vali¬ 
date  the  mathematical  model  were  carried  out  at  full  scale. 

16.  An  instrumentation  system,  designed  to  record  the  motion  of  the  body 
and  described  by  Coman  (Ref  9),  was  fully  contained  within  the  dummy  torpedo 
and  comprised,  three  rate  gyroscopes  to  measure  the  angular  velocity  vector, 
three  accelerometers  to  measure  the  linear  acceleration  vector,  and  a  solid 
state  digital  recorder.  The  trajectory  of  the  body  was  obtained  by  integrating 
the  recorded  angular  velocity  and  linear  acceleration  as  described  by  Coman 
(Ref  10).  The  sensor  signals  were  also  recorded  for  a  ten  second  period 
before  release  and  this  data  was  filtered  to  provide  the  attitude  of  the  body 
at  release  for  the  initial  conditions  of  the  attitude  integration.  The 
initial  conditions  for  the  velocity  integration  were  measured  optically. 

1?.  In  order  to  isolate  the  influence  of  tne  parachute  and  determine  the 
characteristics  of  the  body  alone  the  first  set  of  measurements  were  carried 
out  by  projecting  the  dummy  torpedo  alone, without  any  parachute,  into  the  sea. 

A  second  series  of  measurements  were  then  made  with  parachutes,  the  torpedo 
and  parachute' being  released  from  a  helicopter.  It  was  interesting  to  note 
that  by  commencing  with  the  buoyant  dummy  torpedo,  assumed  stationary,  on  the 
sea  surface  at  the  end  of  a  drop  and  then  integrating  the  measured  motion  back¬ 
wards  through  water  entry  and  through  the  flight  in  air  it  was  possible  t.c 
determine  the  velocity  of  the  delivery  aircraft  and  that  this  value  agreed 
within  1  m/sec  with  the  optically  measured  aircraft  velocity. 

MODEL  VALIDATION 


All 

■-A 

i 


A 


48.  Borne  typical  results  of  the  experimental  measurements  apong  with  the 
predictions  of  the  simulation  are  shown  in  figures  6,  7  and  8.  The  body's 
pitch  angle,  6,  in  degrees  and  axial  component  of  velocity,  u,  in  metres  per 
second  arc  both  plotted  against  time  in  these  figures. 

49.  Figure  6  shows  the  water  entry  behaviour  of  the  bare  torpedo  projected 
into  a  calm  sea  at  approximately  30  m/sec,  at  a  trajectory  angle  of  20°  below 
the  horizontal,  and  with  zero  incidence  to  this  trajectory.  Water  impact 
occurs  at  approximately  0.3  seconds  and  during  the  initial  phases  of  water 
entry  a  nose  down  rate  of  turn  is  imparted  to  the  body  by  the  reduced  pressure 
region  under  the  nose.  At  approximately  0.6  seconds  the  tail  hits  the  top  of 
the  cavity,  in  this  region  the  simulation  diverges  a  little  from  the  measure¬ 
ments  however  this  particular  tail  slapping  behaviour  was  found  to  be  not 
experimentally  repeatable  in  detail. 

60.  In  the  drops  described  in  figures  7  and  8  the  torpedo  was  fitted  with  a 
parachute.  In  figure  7  a  conical  ribbon  parachute  of  approximately  2  metres 
flying  diameter  fitted  with  7.0  metres  long  rigging  lines  was  employed,  the 
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torpedo  was  released  at  a  aeight  of  225  metres  above  sea  level,  water  impact 
occurred  approximately  9*5  seconds  after  release,  and  the  pitch  angle  at  impact 
is  almost  vertical.  In  figure  8  a  ringshot  parachute  of  2  metres  flying  dia¬ 
meter  fitted  with  3.5  metres  long  rigging  lines  was  used,  the  torpedo  was 
released  at  a  height  of  60  metres,  water  impact  occurred  approximately  4  seconds 
after  release  and  the  pitch  angle  at  impact  is  about  55  degrees, 

51.  In  figures  7  and  8  and,  indeed,  in  all  of  the  parachute  drops  which  were 
made  it  was  found  that  the  simulation  predicted  higher  frequencies  of  oscilla¬ 
tion  of  the  body  in  air  then  were  observed.  It  was  not  possible  to  offer  a 
satisfactory  explanation  for  this  inaccuracy  of  the  simulation. 

52.  All  of  the  measurements  which  were  carried  out  in  support  of  the  mathe¬ 
matical  model  described  in  this  note  were  at  impact  velocities  of  between  20 
and  40  m/sec,  however  a  limited  amount  of  work  was  carried  out  to  compare  the 
simulation  with  the  150  m/sec  entry  velocity  full  scale  measurements  described 
by  Waugh  and  Stubstad  (Ref  1,  chap  5).  At  this  higher  impact  velocity  it  is 
to  be  expected  that  the  influence  of  the  underpressure  effect  will  be  reduced. 
Head  shapes  'a',  *g',  * 1  *  and  'n'  were  simulated  and  good  agreement  with  the 
water  entry  whip  (Ref  1,  fig  5-9)  and  with  the  zero  cavitation  number  drag 
coefficient  (Ref  1,  fig  5.1l)  were  obtained  indicating  that  this  simulation  may 
be  applicable  to  a  wide  range  of  impact  velocities. 
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APPENDIX  !5- 

C  FORTRAN  LISTING  OF  WATER  ENTRY  SIMULATION 

C 

C  S.I.  UNITS  ARE  USED 


c 

LENGTH  IN 

METRES 

c 

TIME  IN  SECONDS 

c 

FORCE  IN 

NEWTONS 

c 

MASS  IN  KILOGRAMS 

c 

p 

ANGLE  IN 

RADIANS 

c 

THE  PRINCIPAL  VARIABLES  IN  THIS  PROGRAMME  ARE:- 

c 

NAME  OF  VARIABLE 

DESCRIPTION 

c 

BODY 

PARACHUTE 

c 

VBOD(l) 

VPAR  ( 1 ) 

U  ) 

c 

VBOD  <  2 ) 

VPAR  ( 2  > 

V  )  COMPONENTS  OF  LINEAR  VELOCITY 

c 

VB0D<3) 

VPAR ( 3 ) 

W  ) 

c 

VBOD  (  4 ) 

VPAR(4> 

P  ) 

c 

VBOD <5 ) 

VPAR ( 5 ) 

G  )  COMPONENTS  OF  ANGULAR  VELOCITY 

c 

VBQD(6) 

VPAR (6) 

R  ) 

c 

VBOD <  7 ) 

VPAR ( 7 ) 

EO  ) 

c 

VBOD  ( 8 ) 

VPAR ( 8 ) 

El  )  QUATERNION  ATTITUDE  REPRESENTATION 

c 

VBOD ( 9 ) 

VPAR  ( 9 ) 

E2  ) 

c 

VBOD (10) 

VPAR (10) 

E3  ) 

c 

VBOD (11) 

VPAR (11) 

X  > 

c 

VBOD (12) 

VPAR (12) 

Y  )  POSITION  OF  C  ♦  G «  IN  SPACE 

c 

c 

c 

VBOD < 13) 

VPAR (13) 

Z  > 

XBOD 

XPAR 

) 

c 

YBOD 

YPAR 

)  TRANSFORMATION  MATRIX 

c 

c 

c 

c 

ZBOD 

ZPAR 

) 

ADMBOD 

ADMPAR 

(  THE  FIRST  6  COLUMNS  IS  THE  MASS  MATRIX 
(  THE  7TH  COLUMN  IS  THE  EXTERNAL  FORCES 

c 

p 

BDMBOD 

BDMPAR 

OLD  VALUES  OF  ADMBOD  AND  ADMPAR 

c 

PWBOD  < 1 ) 

PWPAR ( 1 ) 

X/U2  ) 

c 

PWBOD  <  2 ) 

PWPAR ( 2 ) 

Y/V  ) 

c 

PWBOD  <  3 ) 

PWPAR (3) 

Y/R  ) 

c 

PWBOD  <  4 ) 

PWPAR (4 ) 

M/W  )  IN  WATER  FORCE  DERIVATIVES 

c 

PWBOD ( 5 ) 

PWPAR  <  5 ) 

M/G  ) 

c 

PWBOD ( 6 ) 

PWPAR ( 6  > 

X/UDOT) 

c 

PWBOD <  7 ) 

PWPAR ( 7 ) 

Y/VDOT)  #**#*NOTE***** 

c 

PWBOD ( 8 ) 

PWPAR ( 8 ) 

Y/RDOT)  P****(4)=M/W  EXCLUDES  PERFECT 

c 

PWBOD  <  9 ) 

PWPAR ( 9 ) 

N/RDOT)  FLUID  MOMENT  WHICH  IS  CALCULATED 

c 

PWBOD  < 1 0  > 

PWPAR (10) 

Y/V2  )  IN  SUBROUTINE  ENERT 

c 

PWBOD (11) 

PWPAR (11) 

Y/RV  )  I • E ♦ P ( 4 ) =M/W (TOTAL ) +X/UDQT-Y /VDOT 

c 

PWBOD (12) 

PWPAR (12) 

M/W2  ) 

c 

PWBOD (13) 

PWPAR (13) 

M/GW  ) 

c 

PWBOD ( 14) 

PWPAR (14) 

MASS  OF  DISPLACED  FLUID 

c 

PWBOD (15) 

PWPAR (15) 

(DISPLACED  MASS ) # ( X  CO-ORD  OF  C.B.) 

c 

p 

PWBOD ( 16) 

PWPAR (16) 

K/P  ROLL  DAMPING 

L* 

c 

c 

c 

PABOD 

PAPAR 

IN  AIR  EQUIVALENT  TO  PWBOD»PWPAR 

PAPBR 

FULLY  DEPLOYED  VALUES  OF  PAPAR 

c 

PWPBR 

FULLY  DEPLOYED  VALUES  OF  PWPAR 

lb. 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

10 

c 


c 


RBQD 

RPAR 

RELATIVE  MOTION  FIELD 

XAXBOD 

XAXF'AR 

X  CO-ORD 

) 

RADBOD 

RADPAR 

RADIUS 

) 

DAXBOD 

DAXPAR 

DX 

)  EXTERNAL  SHAPE 

SINBOD 

SINF'AR 

SIN( SLOPE ) 

) 

COSBOD 

COSF'AR 

COS (SLOPE) 

) 

HTTBOD 

HTF'AR 

EQUIVALENT 

HEIGHT  OF  ADDED  MASS 

ISW1B0 

)  DEFINES 

TYPE  OF  SEGMENT 

ISW2B0 

) 

ISHD  COUNTS  BROKEN  SHROUD  LINES 

DIMENSION  ADMBOD  <  6  >  7 )  t  BDMBOD  ( 6  » 6 )  » VBQD  ( 13  )  >PABOD<  16)  *PWBOD<  16)  » 
1RB0D (6) » DBOD ( 6) *  XBOD ( 3  > » YBQD ( 3  > » 2B0D( 3 ) » WATER ( 6 ) > 

2XAXB0D( 50 ) i  RADBOD  (  50 ) » DAXBOD ( 50 ) » SINBOD ( 50 ) »C0SB0D(50) » 

3ISW1B0 ( 50 ) f I SW2B0 ( 50 ) »HTTB0D(50) 

DIMENSION  ADMPAR  <  6  r  7  )  t  BDMPAR  ( 6 »  6  >  »  VF'AR  <  13 )  »PAPAR<  16)  »  PWPAR(  16)  i 
1 RPAR ( 6 )  »  DPAR  (  6 ) » XPAR < 3 ) , YPAR ( 3 ) » ZPAR ( 3 ) » PAPBR ( 1 6 ) , PWPBR (16 ) > 
2XAXPARC 10) »RADPAR< 10) rOAXPAR(lO) »SINPAR< 10) t COSPAR < 10)  » 

3 ISHD (30) 

COMMON/BLOCK 1 /  T  t SEAA >  SEAW  t SEAC  » COMX  t COMY  > WDEP  t SURMY » SURMZ » 

1 SGTEM1 » ATSYSZ  >SURECG  >  DURECG  »  BETID »  BET 2D  *  CGOR >  XAXOR  t CAVA  t CAVB  t 
2SGVYVZ» ATVZVY»BETiC»BET2Cf B0FU»B0FV»B0FW»PCAV»XAX»RADrDAX»8INAf 
3CQSA,ISW1 ,  ISW2>  ITES.  ICROS»  IBEP .  IDRY» BODU*  BOBU r  BODW r  PROT 
SY JWEIB.DAT  CONTAINS  INITIAL  CODITIONS 
CALL  ASSIGN( 1 > ' SY ♦ WEIB . DAT '  »  0 » 'OLD' , 'NO' r 1  ) 

8YiUEPB.DAT  CONTAINS  THE  DESCRIPTION  OF  THE  BODY 
CALL  ASSIGN ( 2 » ' SY* WEF'B . DAT ' » 0 » 'OLD' * 'NC' , 1 ) 

SY « OUT  .  DAT  WILL.  CONTAIN  THE  OUTPUT 
CALL  ASS I GN(3> 'SYJOUT.DAT' *0, 'NEW' > 'NC' ,  1  ) 

LENGTH  IN  TIME  OF  SIMULATION 
READ(1»20)TMAX 

TIME  OF  RELEASE  OF  PARACHUTE  < TREL.>TMAXfcNO  PARACHUTE) 

READ  ( 1 »  20 ) TREL 

TIME  OF  COMMENCEMENT  OF  INITIAL  DEPLOYMENT  OF  PARACHUTE 
READ(1 » 20)  TREL 1 

TIME  OP  COMMENCEMENT  OF  DEREEFING  OF  PARACHUTE 
READ ( 1 » 20 ) TREL2 

TORQUE  APPLIED  TO  STORE  BY  RELEASE  LANYARD 
READ ( 1  *  20 ) TORLAN 

TIME  OF  ACTION  OF  TORLAN 
READ  ( 1  >  20 ) TRELAN 

INITIAL  TWIST  IN  SHROUD  LINES 
READ  < 1 »  20 ) TWA 

NUMBER  OF  OUTPUT  RECORDS  REQUIRED 
READ ( 1 »  20 ) PQINUM 
DO 10 1  =  1  *  13 

BODY  INITIAL  CONDITIONS 
READ (1*20) VBOD( I  ) 

WIND  VELOCITY  COMPONENTS  IN  SPACE 
READ (1*20) W INDX 
REAIK1*  20)  WINDY 

SEA  WAVE  AMPLITUDE (METRES) 


gE*^A  -  il  -J  !4*rk  * 


17. 


READ  < 1 » 20 ) SEAA 

C  SEA  WAVE  FREQUENCY (1 /METRES) 

READ  (1*20) SEAW 

C  SEA  WAVE  CELERITY(M./SEC) 

READ  < 1 » 20  >  SEAC 

C  UNIT  VECTOR  IN  DIRECTION  OF  SEAC 

READ  (1*20) COMX 
READ  (1*20) COM Y 

20  FORMAT ( F13 . 5 ) 

21  FORMAT (2F13.5) 

D030I=1 *16 

C  BODY  FORCE  DERIVATIVES 

30  READ  (2*21  )  PABOD<  I )  *  F‘WBQD<  I  > 

C  (TAIL  IN  CAVITY  LIFT  COEFFICIENT) #0.5 

READ (2*20)  CLBOD 

C  (BODY  IN  CAVITY  DRAG  COEFFICIENT ) *0 . 5 

READ (2*20) CDBOD 
C  MASS  OF  BODY 

READ  <  2  *  20  >  PMBOD 
C  Y  CO-ORD  OF  BODY  C.G. 

READ  (2*20)  PY80D 
C  Z  CO-ORD  OF  BODY  C.G. 

READ  (  2 » 20 )  F'ZBOD 

C  BODY  AXIAL  MOMENT  OF  INERTIA 

READ (2*20) PIXBOD 

C  BODY  TRANSVERSE  MOMENT  OF  INERTIA 

READ  (2*20)  F'I  YBOD 
C  MAXIMUM  BODY  RADIUS 

READ (2*20) BODRAD 

C  X  CO-ORD  OF  POINT  OF  ATTACHMENT  OF  PARACHUTE 

READ( 2*20) BODTL 

C  MEASURE  OF  ENERGY  TO  DESTROY  NOSE  CAP 

READ (2» 20) PIMP 
TRAV=0.0 
I PIMP-1 
T=0  ♦  0 

POIDT =TMAX/POINUM 
IFIRST-1 
I SEC0N=1 
D075I-1 *6 
D075J-1 *6 
BDMPARd*  J)=0.0 
75  BDMBQDd*  J>=0.0 

58  NBOD=0 

C  BODY  GEOMETRY  TEM2=X  CO-ORD*  TEM4=RADIUS* 

C  TEM8=ADDED  MASS  HEIGHT 

C  11=0  AND  Jl=l  ORDINARY  BUOYANT  SECTION 

C  11=1  AND  Jl=l  FLOODED  SECTION  EG  PARACHUTE  CONTAINER 

C  11=1  AND  J1=0  CRUCIFORM  TAIL 

C  11=0  AND  J1=0  SHROUD  RING  TAIL 

C  SET  TEM4=-1 .0  AT  END  OK  COMPLETE  BODY  DATA  BLOCK 

C  TWO  DATA  BLOCKS  ARE  REQUIRED  FIRST  WITH  NOSE  CAP 

C  SECOND  WITHOUT  NOSE  CAP 

C  IF  PIMPOO.O  THE  FIRST  BLOCK  IS  IRRELEVANT 

60  READ (2*62) TEM2  * TEM4 * TEM8* 1 1  *  J1 


xo. 


42  FORMAT (3F13«5»2I2) 

IF ( TEM4 .  LT  *  -0 . 1 )G0T02O0 

IF  < 1 1 . NE « I ♦ OR . J1 . NE ♦ J ) 60T0A7 

TEM5=TEM1-TEM2 

TEM6=TEM4-TEM3 

TEM7=SGRT ( TEM5*TE MSTTEM4*TEM6 ) 

IF<  <  TEM5/TEM7 ) . GT . 0 ♦ 002 ) G0T065 
TEM5=0 , 002*ABS<  TEM6 > 

TEM7=SQRT ( TEM5*TEM5+TEM6*TEM6 ) 

45  NBQD^NBQD+l 

XAXBOD <  NBOB )  =  <  TEM1+TEM2 > /2 . 0 
RABBOB  <  NBOD )  =  f TEM3+TEM4 ) /2  4  0 
HTTBOD ( NBOB ) = ( TEM8+TEM9 ) /2 ♦ 0 
DAXBOB ( NBOB ) =TEM5 
SINBOD<  NBOB ) =TEM4/TEM7 
COSBOB ( NBOB )=TEM5/TEM7 
ISW1B0 ( NBOB ) = I 1 
ISU2B0  <  NBOB ) = J1 
47  1=11 

J=J1 

TEM1=TEM2 

T£M3=TEM4 

TEM9-TEM8 

G0T040 

4000  ISEC0N=0 

C  PARACHUTE  INITIAL  CONBITIONS 

VPAR<1  )«VBOD<  1  > 

VPAR  (  2  )  =VBOD  <  2  )  +BODTL.*VBQD  (  4  ) 

VPAR ( 3 ) =VBOD  <  3 ) -BODTL*VBOD < 5 ) 

DQ6050I=4 » 10 
4050  VPAR ( I ) =VBQD  < I > 

VPAR <11 >*VB0D<11)+B0DTL#XBGD<1> 

VPAR < 12 ) =VBGD (12) +BOBTL*XBOB <  2 ) 

VPAR  < 13 ) =VBOB <13) +BODTL*XBOD  <  3  > 

C  SY ♦ WEPP ♦ BAT  CONTAINS  THE  DESCRIPTION  Or  THE  PARACHUTE 

CALL  ASSIGN  <  4  * ' SY l WEPP «  BAT ' rOr ' OLB ' > 'NC' » 1 ) 

C  NUMBER  OF  SHROUB  LINES 

READ  <  4  >  50 ) NSHB 
V50  FORMAT  <  12) 

C  TORSIONAL  STIFFNESS  OF  SHROUB  LINE  SYSTEM 

REAB ( 4  »  20 ) TUI  ST 

C  LENGTH  OF  EACH  SHROUB  LINE 

READ (4 » 20 ) TLFL 

C  TEAR  STRIP  YIELB  LOAD 

READ  (  4  #  20  )  TL.FT 

C  TEAR  STRIP  EXPIRED  LENGTH  +  TLFL 

REAB  <4*20) TLFTL 

C  TEAR  STRIP  BREAKING  LOAD 

REAB (4  t 20) TLFF 

C  TLFA  f  TL.FB » TL.FC t  TLFP  DESCRIBE  SHROUB  LINE  STRESS/STRAIN 

C  CHARACTERISTIC  SEE  7000 

REAB  <  4  >  20 ) TLFA 
REAB ( 4  »  20 ) TLFB 
REAB  <  4  »  20 ) TLFC 
READ <  4  >  20 ) TLFP 


C  IN  WATER  VALUE  Of  TLFC 

READ<  4f 20)TLCWET 
C  IN  WATER  VALUE  OF  TLFP 

READ  (4*20)  TLF'WET 

C  SHROUD  LINE  BREAKING  LOAD 

READ (4*20) TLFBR 

C  RATE  OF  DEPLOYMENT  PARAMETER 

READ  (4*20) SDLPAR 
D02030I=1 *16 

C  PARACHUTE  FORCE  DERIVATIVES 

2030  READ  <  4 , 21 ) PAPBR (I )  * PWPBR  < I ) 

C  PARACHUTE  ADDED  MASS  HEIGHT 

READ  (4*20)  HTPAR 

C  PARACHUTE  PARTIALLY  WET  DRAG  CQEFFICIENT*0 . 5 

READ  (4*20) CDPAR 
C  MASS  OF  PARACHUTE 

READ <4  *  20) PMPAR 

C  PARACHUTE  AXIAL  MOMENT  OF  INERTIA 

READ (4*20)  PIXPAR 

C  PARACHUTE  TRANSVERSE  MOMENT  OF  INERTIA 

READ  (4*20) PI YPAR 

C  MAXIMUM  RADIUS  OF  FULLY  DEPLOYED  PARACHUTE 

READ  <  4  »  20 )  F'BRRAD 

C  X  CO-ORD  OF  POINT  OF  ATTACHMENT  OF  SHROUD  LINES  TO  CANOPY 

READ (4*20) PARTL 

C  RADIUS  OF  PARACHUTE  AT  RELEASE  AS  FRACTION  OF  PBRRAD 

READ (4*  20) SCLMIN 

C  RADIUS  OF  DROGUE  AS  FRACTION  OF  PBRRAD 

READ ( 4 » 20 ) SCLMAX 

C  RADIUS  OF  REEFED  PARACHUTE  AS  FRACTION  OF  PBRRAD 

READ <4 *  20 ) SCL2 
NPAR»0 

C  FULLY  DEPLOYED  PARACHUTE  GEOMETRY  TEM1=X  CO-ORD *  TEM3=RADIUS 

READ <4 *  2002 ) TEM1 *  TEM3 
2000  RE  ADM*  2002 )TEM2*TEM4 
2002  FORMAT  <  2F13 « 5 ) 

IF  <  TEM4 . LT . ~0 ♦ 1 ) G0T02060 
TEM5=TEM1 -TEM2 
TEM6=TEM4-TEM3 

TEM7=SQRT <  TEM5*TEM5+T£M6*TEM6 ) 

IF < (TEM5/TEM7) .GT. 0.002 )G0T02005 
TEM5=0.002*ABS(TEM6) 

TEM7=SQRT ( TEM5*TEM5+TEM6#TEM6 ) 

2005  NPAR=NPAR+1 

XAXPAR ( NPAR )  =  <  TEM 1+TEM2 ) 72 ♦ 0 

RADPAR ( NPAR  >  =  ( TEM3+TEM4 ) /2 . 0 

DAXPAR ( NPAR ) =TEM5 

SINPAR  <  NPAR  >=TEM6/T£M7 

COSPAR < NPAR >=TEM5/TEM7 

TEM1=TEM2 

TEM3=T£M4 

G0T02000 

2060  SHDANG=6 . 283185/FLOAT  <  NSHD ) 

D02065I-1 *NSHD 
2065  ISHD ( I ) =0 


20. 


C 

c 

200 

202 

205 

C 


C 


C 


C 


208 

C 

210 


C 


SCLPAR=SCLMIN 

dt*o.o 

GOT06005 

FROM  200  TO  1500  IS  FORCES  ON  BODY 
TEST  FOR  PRESENCE  OF  NOSE  CAP 
IF ( IPIMP «  EG «  0 ) GOT0202 
IF ( TRAV*VEL2 . LT . PIMP ) G0T0202 
IP IMP-0 
G0T058 
D0205I=1»6 
D0205J=1 » 7 
ADMPAR  < I » J ) =0 , 0 
ADMBOB ( I » J ) =0 . 0 
CALL  XY2 ( VBOD *  XBOD »  YBOD »  ZBOD  > 

CALL  SEA  <  WATER  >  XBOD  t YBOD  *  ZBOD » VBOD  <  1 1 )  , VBOD  < 1 2 ) , VBOD (13)) 

DEPTH  OF  EXTREMITIES  OF  BODY 
DEP=SURECG-XAXBOD ( 1 ) *DURECG 
TEM2=SURECG-XAXB0D ( NBOD ) *DURECG 
DEPME=  ( DEP+TEM2 )  /2  4  0 
TEM7=B0DRAD*SGTEM1 
TEM3=DEP-TEM7 
TEM4=T£M2-TEM7 
TEM1=DEP+TEM7+B0DRAD 
TEM2=TEM2+TEM7+B0DRAD 

FIRST  ESTIMATE  OF  TIME  INCREMENT 
DT-BODRAD/SQRT < VBOD ( 1 ) *VBGD ( 1 ) +VBQD<  2 ) #VBOD ( 2 ) +VBOD ( 3 ) *VBQD ( 3 ) ) 

I F  <  DT ♦ GT  *  0 . 02  > DT=0 . 02 
I F ( DEPME . GT . 0 . 0 ) G0T0208 

RELATIVE  MOTION  FIELD  IN  AIR 
RBODU  >“VB0D<1 )-WINDX*X»OD< 1 >-WlNDY*XB0D<2> 

RBOD < 2 ) =VBOH <  2 ) -WINDX*YBOD < 1 ) -WINDY*YBOD < 2 > 

RBOD ( 3 ) =VB0D ( 3 ) -W INDXSZBOD ( 1 ) -WINDY*ZBOD ( 2  > 

RBOD  <  4 ) =-9  4  81*XBGD <  3 ) 

RBOD(5)=-9.31*YBOD<3) 

RBOD ( 6 ) =-9. 81#ZB0D<  3 ) 

CALL  HYDRO ( ADMDOD > PABOD » VBOD » RBOD ) 

I F ( T , LT . TRELAN ) ADMBOD ( 5 » 7 ) =ADMBQD ( 5 1 7 ) +TORLAN 
I F ( TFM1 , GT 4  0 ♦ 0, OR . TEM2 4 GT 4 0 . 0 )60T0210 
BODY  IS  FULLY  IN  AIR 
PCAV=101000  4  0 
TRAV=0 . 0 
1 1 M = 0 
GOT0400 
1 1  M=  1 

RELATIVE  MOTION  FIELD  IN  WATER 

RBOD ( 1 ) =VBOD ( 1 ) -WATER  < 1 ) *XBOD ( 1 ) -WATER ( 2 ) *XBOD ( 2 ) -WATER ( 3 ) *XBOD ( 3  > 
RBOD  < 2 ) »VBQD <  2  > -WATER  < 1 ) A YBOD ( 1 ) -WATER ( 2 ) SYBOD ( 2 ) -WATER ( 3 ) # YBOD ( 3 ) 
RBOD  ( 3 )  =VBOD(  3 )  -WATER <  1  >  #ZBQD  ( 1 ) -WATER  (  2 )  #ZB  JiH  2  > -WATER  ( 3 )  JKZBOD  ( 3 ) 
RBOD ( 4 ) =WATER  <  4 ) *XBOD ( 1 )+ WATER <  5 ) #XBOD ( 2 )  +  <  WATER < 6 ) -9 . 81 ) *XBOD ( 3 > 
RBOD<5>=WATER<4)*YBOD<1)+WATER<5)*YBOD<2)+(WATER<6>-9.31 >*YB0D<3) 
RB0D(6>=WATER<4)*ZBQD<1 >+WATER(5)*ZB0D(2)  «•<  WATER (6) -9,81 )*ZB0D(3) 
OBTAIN  CAVITATION  NUMBER 
V2W2-RB0D ( 2 ) *RBQD ( 2 ) +RBOD ( 3 ) #RBGB <  3 ) 

VEL. 2=V2W2+RB0D(  I  >*RBOD<  1 ) 

PAMB=10100. 0*WDEP 


r-o  n 


IF(PCAV4EG.0.0)G0T0216 
PCAV=PAMB-TRAV*300  4  O/BOBRAD 
IF  <  PCAV  *  LT . 0 . 0 ) PCAV«0 . 0 
PAMB=PAMB-PGAV 
216  CAV=PAMB/<515,0*VEL2> 

I F  <  TEM3 ♦ LT ♦ 0  4  0  4  OR . TEM4 . LT  4  0 . 0 ) G0T0240 
C  BODY  IS  FULLY  IMMERSED  IN  WATER 

IF(CAV,LT4043)G0T0250 
C  THERE  IS  NO  CAVITY 

CALL  HYDRO  <  ADMBOD » PWBOD » VBOD  t  RBOD ) 

TRAV=9999 , 0 
PCAV=0 . 0 
GOT0400 

BODY  IS  PARTIALY  WET 
40  DT=DT/5.0 

250  ICR0S=1 
ICS=0 
IUP«0 

SINMAX=SGRT < V2W2/VEL2 ) 

IF  <  DEPME . GT . 0 . 0 ) CALL  HYDRO  <  ADMBOD » PABOD »  VBOD 4  RBOD ) 

I F ( 1 1 M ♦ EG ♦ 0 . OR . PCAV , EG . 0 . 0  >  G0T0256 
C  TRAV=D I STANCE  TRAVELLED  AFTER  IMPACT 

TEM8-VB0D (11) -TORX 
T£M9=VB0D( 12) -TORY 
TEM10=VB0D(13>-T0RZ 

TRAV=SGRT<TEMQ*TEM8+TEM9*TEM9+TEM10#TEM10> 

C  CALCULATE  UNDERPRESSURE 

256  DUNF'E*  101000  4  0*  <  1 . 0-  ( 0 . 2-0 . 07*XB0D  ( 3  >  >  #TRA V/BODRAD ) 

C  CALCULATE  FORCES  AND  ADDED  MASSES  FOR  EACH  BODY  SEGMENT 

DQ1500J=1 f NBOD 
I TES-0 

XAX=XAXBOD< J) 

RAD«RADBOD< J) 

DAX=DAXBOD<  J) 

SINA=SINBOD(J) 

C0SA=CQSBQD  <  J ) 

ISW1=ISW1B0< J) 

ISU2=ISW2B0( J) 

PROT =VBQD (4 ) 

BODU“RBQD  < 1 ) 

B0DV=RB0D(2)+XAX*VB0D<6> 

B0DW  =  RB0D <  3 ) -XAX*VBOP <  5 ) 

IF ( I UP . EG . 0 ♦ OR 4  X AX . GT . CGOUP )G0T01110 
C  UPWARD  VELOCITY  ON  REAR  OF  CAVITY 

TEM1 = ( CGQUP-XAX ) /CAVUP 
IF<TEM1 .GT. 1,0>TEM1-1.0 
B0DV=B0DV+TEM1*UPVEL 
B0DW=B0DW+TEM1*UPWEL 
1110  B0FU=B0DU4SINA 

B0FV=BGDV#CQSA-RAD#VB0D<6)#SINA 
&0FW=B0DW*C0SA+RAIi*VBQD(5>#SINA 
CALL  DEPTH 

IF< I  DRY. EG. 1) GOTO 1500 
IFdIM.EQ.i ) GOTO 1 120 
SET  UP  IMPACT  POSITION 


C 


I  IM=1 

T0RX=VB0D<11> 

TORY=VBQD (12) 

T0RZ=VB0D< 13) 

C  IF  ICR0S=1  AXIAL  ELLIPSOID  CAVITY  NOT  FORMED 

C  IF  ICR0S*=0  AXIAL  ELLIPSOID  CAVITY  EXISTS 

.1120  IF<ICROS.EQ,0>GOT01170 
IFaSW2.EG.0>G0T01150 
C  TEST  SEGMENT  FOR  CAVITATION 

C  BET 1C  TO  BET2C  WILL  BE  THE  WETTED  SECTOR  OF  THE  SEGMENT 

1130  TEM7=B0DV#B0DV4BQDW*B0DW 
U2V2W2=TEM74B0DU*B0DU 
IF<F'CAV ♦  GT ♦ 0 »  0 )  GOTO  1 132 
PAMB=  10 1 00 , 0*  ( WDEP-XAX*RBOD  ( 4  ) /? .  81 ) 

1132  CAVSEG=F'AMB/  (515, 0*U2V2W2 ) 

IF<  CAVSEG • GT , 1 ♦ 0 )G0T01 150 

IF ( CAVSEG » LT , 0 , 001 ) CAVSEG=0  *  001 

TEM3=SGRT (U2V2W2) 

TEM27=SGRT(TEM7> 

TEM1=CGSA*TEM27 

TEM4=<B0FU4TEM1)/TEM6 

IF'  ( SINMAX  » LT ,  TEM4 ) SINMAX=T£M4 

J1==J-1 

IF(J.EG»1)J1-1 

C  CALCULATE  INCIDENCE  CONDITION  FOR  CAVITATION 

TEM28=RAD*<SINB0D< J1 > *C0SA-C0SB0D< J1 >*SINA> 

TEM28=TEM28* ( COSBOD ( J 1 ) +CQSA ) / <  DAXBOD ( J1 ) 4DAX ) 

T£M2=0,2S*TRAV/RAD 

I F ( TEM2 . GT . 0 . 76 ) TEM2=0 . 76 

T£M2~TEM2*(0.240.S*SINMAX*SINMAX> 

CaNlNC=TEM2#SlNMAX-0,34*EXP<-TEM28>-0,36*CAVSEG4Q,28*TEM7/U2V2W2 
TEM2~CQNINC*TEM4 
IF< (B0FU4TEM1) . GT . TEM2 ) G0T01 1 40 
C  SEGMENT  IS  ALL  DRY 

C  CALCULATE  GEOMETRY  OF  AXIAL  CAVITY  IF  PRESENT 

T£M2=TEM27/TEM6 
TEM3-ABS < BODU ) /TEM6 
TEM4=TEM3-ABS(TEM2#SINA/C0SA) 

IF'  (  TEM4 .  LE  ,  0 , 0  >G0T01400 
TEM1=RAD/TEM4 

CAVA=0,4*SINMAX*<1 , 04TEM28 ) #TEM1 /CAVSEG 
CAVB=TEM140, 13#CAVA 
TEM9-1  s'5*TEM140.09#TRAV 
IF ( TEM9  *  LT ♦ CAVB)CAVB=TEMV 
ICS  =  0 

l'F  (  (  CAVA4CAVA )  ,  LT .  TRAV ) GOTO  1 1 35 

CAVA~0,84*TRAV 

CAVC=0.21#TRAV 

ICS=1 

1135  IF (CAVA,L£,0,0)GOT01400 
C  FIT  CAVITY  TO  NOSE 

TLM4-TEM1 /CAVB 

CX0R=CAVA*SQRT  ( 1 . 0-TEM4*T£M4  ) 

TEM4=-TEM4#CAVA*CAVA/CAVB 

TEM5=CX0R/SGRT(CX0R*CX0R4TEM4*TEM4) 


23. 


I 


i 


TEM4=TEM2*C0SA+ABS  <  TEM3KSINA ) 

IF ( TEM5 .LE  4  TEM4 ) GOTO 14 00 

CAVY=-B0DV/TEh6 

CAVZ*-B0DW/TEM6 

CGOR=CXOR-XAX 

XAXOR=XAX 

SGVYVZ=SGRT  <  CAVY*CAVY+CAVZ*CAVZ ) 

ATVZVY«=ATAN2  <  CAVZ  1  CAVY ) 

ICR0S*0 

IF  <  XAX .LT«0*0«0R«ICS.EG»1> GOTO 1400 
C  UPWARD  VELOCITY  ON  REAR  OF  CAVITY 

CGOUP=-CGOR 
CAVUP=CAVA 
IUP»1 

UPVEL=-TEM6* < 0 . 005*RBOD ( 5 ) +31 GN ( 0 ♦ 025 » RBQD ( 5  > )  ) 
UPWEL=-TEM6# ( 0 . 005*RBQD <  6 ) +SIGN  <0 ♦ 025  t RBOD  <  6 ) > ) 
IF<ABS<XB0D<3 ) ) . LT . 0 . 96>G0T01400 
UPVEL  =  TEM4*SIGN<0.025>VBQD<6>  > 
UPWEL=-TEM6*SIGN<0.025f VB0D<5) ) 

G0T01400 

1140  IF  < <  B0FU-TEM1 ) ♦ GE . TEM2 ) G0T01 150 

C  SEGMENT  IS  PARTIALLY  WET  BY  CAVITATION  CONDITION 

TEM2= <  TEM2-B0FU ) /TEM1 
TEM1=ATAN2<  BODV » BODW ) 

TEM3=SQRT< 1 . 0~TEM2*TEM2 > 

TEM4=AT  AN2  <TEM2»TEM3) 

BET  1C=TEM4--TEM1 
BET2C=3, 1416-TEM4-TEM1 
G0T01200 

C  SEGMENT  IS  COMPLETELY  WET  BY  CAVITATION  CONDITION 

1150  BET1C=0 ♦ 0 

BET2C=6  *  2831 

ITES=1 

GOT01200 

C  TEST  FOR  INTERSECTION  WITH  CAVITY 

1170  TEM1=CG0R+XAX 
TEM8-CAVA 

IF  <  TEM1 . GT  * 0 . 0 , OR . ICS,EG.0>GQT01175 
TEM8=CAVC 

1175  TEM1=TEM1/TEM8 
TEM1=TEM1*TEM1 
IF  <  TEM1 . GE . 1  *  0 ) G0TQ1 1 80 
CAVR=CAVB*SQRT< 1.0-TEM1 ) 

CAVD-ABS <  XAXOR-XAX ) *SGVYVZ 

TEM1-CAVR+RAD 

IF  <  CAVD . GE . TEM1 ) GOTO 1 180 

TEM1-CAVR-RAD 

IF ( CAVD . LE ♦ TEM1 ) GOTO 1400 

TEM1 =-TEMl 

I F  < CAVD . LE »  TEM1 ) G0T01 180 

TEM1  =  <  CAVD*CAVD+RAD*RAD-CAVR*CAVR ) / <  2 . 0*CAVD*RAD ) 
TEM2-SQRT  <1.0-TEM1#TEM1 ) 

TEM3=ATAN2<TEM2»TEM1 ) 

BET1C=ATVZVY+T£M3 

BET2C=ATVZVY--TEM3 


l, 


LA 
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:L  180 


1200 

1400 

C 


1500 

400 

C 

C 


C 


c 

c 

3  2 1 0 


C 

c 


GQTQ1200 

IF<ISU2,EG,0>GOTQ1150 

ICR0S“1 

S INMAX“SGRT ( V2W2/VEL2 ) 

G0T01130 

CALL  FORCES ( ADMBOD » RBOD t HTTBQD ( J ) » CLBOB » COBOD ) 

IF< ICRQS.EG. 1 .AND. ICS . EG  *  0) GOT 01500 
IF  < IDER . EG . 1 . OR .  DUNPE  , LT . 0 ♦ 0 . OR ,  ISW2  ,EQ,0)  G0T01 500 
UNDERPRESSURE  FORCE 
TEM4=DUNPE*DAX*RAD 
TEM1=TEM4* (COS (BETID) -COS < BET2D  >  ) 

TEM2=TEM4* < SIN (BET2D) -SIN C BETID)  ) 

TEM3-XAX~RAD*SINA/CQSA 
ADMBOD (2,7) “ADMBOD (2,7 ) +TEM2 
ADMBOD (3,7) “ADMBOD ( 3 , 7 ) +TEM 1 
ADMB0D<5f 7)=ADMB0D<5,7)-TEM3*TEM1 
ADMBOD (6,7) “ADMBOD (6,7 ) +TEM3*TEM2 
CONTINUE 

IF ( ISECON « EG • 1 > G0T05500 

UP  TO  4500  IS  FORCES  ON  PARACHUTE 
CALL  XYZ ( VPAR  ,  XPAR , YPAR »  ZPAR ) 

CALL  SEA ( WATER , XPAR  » YPAR , ZPAR , UPAR ( 1 1  )  »  VPAR ( 12 )  >  VPAR ( 13 ) ) 

DEPTH  OF  EXTREMITIES  OF  PARACHUTE 
0£P= SURE CG-XA XPAR ( 1 ) *DURECG 
T£M2“SURECG-XAXPAR<NPAR)*DURECG 
TEM7“F'ARRAD*SGT£Ml 
TEM1  =  DEF'+TEM7 
TEM13»DEP-TEh? 

TEM14-TEM2-TEM7 
DEPME“  <  TEM1+TEM14 ) /2. 0 
IF ( DEPME , GT . 0 . 0 ) G0T0321 0 

RELATIVE  MOTION  FIELD  IN  AIR 
RPARC 1 >-VPAR< l)-WINDX*XPAR<l)-WINDY*XPAR<2> 

RPAR ( 2 ) * VPAR ( 2 ) -WINDX* YPAR ( 1 > -WINDY* YPAR ( 2 ) 

RPAR <  3 ) -VPAR ( 3 ) -WINDX*ZPAR (1)-UINDY*ZPAR(2) 

RPAR ( 4 ) =~9 .  81 *XPAR ( 3 ) 

RPAR  <  5 )  --■? .  8 1  *  YPAR  ( 3 ) 

RPAR<A)=-9.81*ZPAR<3) 

call  HYDRO  <  ADMPAR > PAF'AR  >  VPAR  >  RPAR  ) 

IF(TEMl.LT.O.O) G0TO3400 
PARACHUTE  COULD  BE  WET 
RELATIVE  MOTION  FIELD  IN  WATER 

RPAR  < 1 >*VPAK< 1 ) “WATER ( 1 )*XPAR( 1 ) -WATER ( 2 ) *XPAR ( 2 ) -WATER ( 3 ) *XPAR ( 3 ) 
RPAR < 2 ) = VPAR v  2 ) -WATER ( 1 ) * YPAR < 1 > -WATER ( 2 ) * YPAR ( 2 ) -WC . ER < 3 ) *YPAR ( 3 ) 
RPAR  (  3 )  -  VPAR  (  3 )  -WATER  ( 1 )  *ZPAR  ( 1 )  -  WATER  <  2 )  *ZF’AR  ( 2 )  -WATER  ( 3 )  *ZPAR  <  3 ) 
RPAR ( 4 ) “WATER ( A ) *XPAR < 1 )+ WATER ( 5 ) *XPAK ( 2 )  +  ( WATER  <  6 ) -9 . 8 1 ) *XPAR ( 3 ) 
RPAR  <  5 )  =  WATER  <  4  )  *YPAR  <  1 )+ WATER  <  5 )  *YPAR <  2  >  +  <  WATER  (  6 )  -9 . 81  >  *YF'AR  ( 3 ) 
RPAR  C  6 )  -WATER  ( 4 )  *ZF'AR  <  1 )  +WATER  (  5  )  #ZPAR  (  2 )  +  <  WATER  <  6  >  ,  81  >  *ZF'AR  (  3 ) 

IF < TEM13 . GT , 0,0.  AND , TEM1 4 , GT . 0 , 0 . AND . CAV . GT , 0 ♦ 3 ) G0T03240 
PARACHUTE  COULD  BE  PART  WET 

SET  UP  TEST  FOR  INTERSECTION  WITH  CAVITY  FROM  BODY 
I F ( I CRQS »  EG . 0 ) G0T03220 
IF  <  CAV. LT ,0.001 )CAV  =  0,001 
CAVA=2.0#B0DRAD/CAV 
CAVB=B0DRAD+0 . 13*CAVA 


3220  TEM1-T0RX-UB0D (11) 

TEM2®T0RY-UB0D<12) 

TEM3=T0RZ~0B0D<13) 

TEM15=TEM1#XF'AR<  1  )+TEM2*XPAR<2)+T£M3#XPAR<3) 

IF ( TEM15. EQ . 0 . 0 ) G0T03230 
TEM4«UB0D < 1 1 ) -  VPAR (11) 

TEM5»VB0D  < 1 2 ) -UPAR  < 1 2 ) 

TEMA=UBOD  < 13 ) -VPAR ( 13 ) 

TEM16-TEM4*XFAR ( 1 ) +TEM5*XFAR ( 2 ) +TEM6KXPAR  <  3  ) 

TEM17=-TEM16/TEM1 5 

IF  <  TEMI  7  4  LE  *  0  *  0 )G0T03230 

7EM18«SGRT<TEM1*TEM1+TEM2*TEM2+TEM3*TEM3) 
7EM1=T£M4+TEM1#T£.  1? 

T£M2=TEM5+TEM2*TEM17 

TEM3*T£M6+TEM3#TEM17 

CAUY=TEM 1  *YPAR  ( 1 )  +TEM2*  YPAR  ( 2 )  +TEM3*YPAR  ('  3 ) 
CAUZ=T£M1*ZPAR<1 > +TEM2*ZPAR < 2 ) +TEM3*ZPAR < 3 ) 
CAVD=SGRT<CAUY*CAUY-rCAV7*CA0Z) 

TEM1=1 .0-TEM17#TEi118/CAVA 
TEM1=TEM1*TEM1 
IF ( TEMI . GE « 1 . 0 ) GGT03230 
CAUR=CAUB#SQRT  < 1 . 0-TEK1 > 

TEM1=CAVR+PARRAD 

IF<  CAVD . GE . TEM1 ) GOT03230 

TEM1-CAVR-PARRAD 

IF  <  CAVD « LE  ♦ TEMI ) G0T03225 

ATUZVY=ATAN2<CAVZ  fC  'Y) 

ICAV=1 

G0T03250 

3225  IF  ( DEPME .  GT « 0 . 0  )  CALL  HYDRO  ( ADMPAR f PAPAR  f  VPAR  f  RPAR ) 
G0T03400 
3230  ICAV=0 

IF ( TEMI 3 , LT . 0 . 0 . OR . TEM 1 4 . LT . 0 . 0 )  GOT03250 
C  PARACHUTE  IS  FULLY  WET 

3240  CALL  HYDRO < ADMPAR f PWPARf VPAR f RPAR) 

G0T03400 

3250  IF  < DEPME . GT , 0 . 0 ) CALL  HYDRO ( ADMPAR .. PAPAR  f  VPAR  f  RPAR  ) 

C  CALCULATE  FORCES  AND  ADDED  MASSES  ON  EACH  SEGMENT 

D04500  J=  1  f  NF'AR 
ITES=0 

XAX-XAXF'AR  <  J ) 

RAD=SCLPAR*RADPAR <  J ) 

DAX“SCLF’AR*DAXPAR  <  J ) 

SINA-SINF'AR  ( J ) 

COSA=COSPAR< J) 

ISW1=1 

ISW2=1 

B0DU-RPAR< 1 > 

D0FU=BQDU*SINA 
TEM1=XAX*C0SA-RAD#SINA 
LOFV=RPAR  <  2  )  *COSA+VPAR  (  6  )  *TEM  1 
DOFW=RPAR ( 3 ) #COSA- VPAR ( 5 ) #TEM  1 
CALL  DEPTH 

IF< IDRY . EG . 1 )G0T04500 
IF i ICAV. EG. 0>GCT04350 
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TEM1=CAVR+RAD 
IF(CAVD*GE*TEM1  )G0T04350 
T£M1«CAVR-RAD 
IF(CAVD,LE,TEM1 >G0TQ4500 
TEM1=-TEM1 

IF (CAVD«LE«  TEH1 )G0T04350 

TEMl=a  ( CAVD*CAVD+RAD*RAD-CAVR*CAVR )  /  <  2 . 0*CAVD*RAD ) 
TEM2=SQRT<1 .0-TEMl*TEMl > 

T£M3*=ATAN2(  TEM2»  TEM1 ) 

B£T1C=ATVZVY+T£M3 

BET2C=ATVZVY~TEM3 

GOTG4400 

C  ALL  WET  BY  CAVITY  INTERSECTION 

4350  BET1C=0.0 

BET2C=6  *  2831 
ITES^l 

4400  CALL  FORCES  <  ADMPAR  t  RPAR » HTF'AR »  0 . 0 » CDF’AR ) 

TLFC-TLCWET 
TLFP=TLPWET 
4500  CONTINUE 

C  CALCULATE  FORCES  IN  SHROUD  LINES 

C  FIRST  OBTAIN  BODY  TO  PARACHUTE  TRANSFORMATION  MATRIX 

3400  S1=XB0H<  1  )*XPAR(1  > +XB0D ( 2 ) *XPAR ( 2 ) +XBOD < 3 ) *XPAR ( 3 > 

S2=YB0D  (  1 )  KcX PAR  ( 1 )  +YBOD  <  2 )  *XPAR  ( 2 )  +YBOD  ( 3 )  *XPAR  ( 3  > 

S3=ZB0D ( 1 ) *XPAR ( 1 ) +ZBOD  <  2  >  #XPAR  <  2 )  +  ZBOD  <  3 ) *XPAR ( 3 ) 

S4=XB0D < 1 ) #YPAR < 1 ) +XBOD ( 2 ) *YPAR  <  2 ) +XBOD <  3 ) *YPAR ( 3  ) 
S5»YB0D(1 )*YPAR( 1 > +YB0D ( 2 > *YPAR ( 2 > +YBQD ( 3 > *YPAR ( 3 ) 
S6=ZB0D(1 )*YPAR(1 ) +ZB0D ( 2 ) *YPAR ( 2 > +ZB0D ( 3 > *YPAR < 3 > 

S7=XB0D ( 1 ) *ZPAR ( 1 ) +XBOD ( 2 ) *ZPAR ( 2 ) +XBOD C  3 ) *ZPAR ( 3 ) 

S8-YB0D ( 1 ) *ZPAR ( 1 ) +YBOD <  2 > *ZPAR ( 2 ) +YBOD ( 3 ) *ZPAR  <  3  ) 
S9=>ZB0D< 1 >*ZPAR( 1 ) +ZBQD ( 2 > *ZPAR ( 2 > +ZB0D ( 3 > *ZPAR < 3 > 

C  RELATIVE  VELOCITY  AND  POSITION  OF  BODY  W.R.T.  PARACHUTE 

T EMI -VBOD ( 2 ) +VBOD ( 6 ) #BODTL 
TEM2= VBOD <  3 ) -VBOD ( 5 ) *BODTL 
B0DU«*VB0D ( 1 ) *S1 +TEM1 *S2+TEM2*S3 
BODV=VBOIi(l  >*S4  +  TEM1*S5+TEM2*S6 
BOD W~ VBOD ( 1 ) #S7+TEM1 #S8+TEM2#S9 
TEM1 »VBOD (11) -VPAR  (11) 

TEM2=VBQD  < 1 2 ) -VPAR  < 12 ) 

TEM3=VB0D (13) -VPAR (13) 

TLX=TEM1*XPAR< 1 > +TEM2*XPAR < 2 > +TEM3*XPAR ( 3 ) +B0DTL*S1 
TLY=TEM1*YPAR< 1 > +TEM2*YPAR ( 2 ) +TEM3*YPAR ( 3 ) +B0DTL*S4 
TLZ-TEM1*ZPAR(1 > +TEM2*ZPAR ( 2 > +TEM3*ZPAR < 3 > +BODTL*S? 

C  CALCULATE  AND  SUM  TENSIONS  IN  SHROUD  LINES 

7000  TPX-0.0 

TPY=0  *  0 
TPZ-0 ♦ 0 
TEM12--SHDANG 
D07500I=1 t NSHD 
TEM12=TEM12+SHDANG 
IF( ISHD( I ) .EG. 1 >GQT0?450 
C  EXTENSION  OF  LINE 

TEM1 -TLX-PARTL 
TEM6=PARRAD*0QS ( TEM12 ) 

TEM2-TL.Y-TEM6 
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TEM7=PARRAD#SIN(TEM12) 

TEM3-TLZ-TEM7 

TEM4=SGRT(TEM1*TEM1+TEM2*TEM2+TEM3*TEM3> 
TEM5-TEM4-TLFL 
IF ( TEM5 . L£ . 0 . 0 ) 60T07500 
C  RATE  OF  EXTENSION  OF  LINE 

TEM8-B0DU- <  VPAR  < 1 ) +VPAR ( 5 ) *TEM7- VPAR  <  6 ) *TEM6 ) 
TEM9-B0DV-  ( VPAR  ( 2  >  +VPAR  ( 6  >  *PARTL-VPAR  <  4  )*TEM7 ) 
TEM10=B0DW-  ( VPAR  ( 3 )  +VPAR  ( 4 )  #TEM6- VPAR  (  5 ) *PARTL ) 
TEM1-TEM1/TEM4 
TEM2-TEM2/TEM4 
TEM3=TEM3/TEM4 
C  TENSION  IN  LINE 

TEM4=(TEM8*TEM1+TEM9*TEM2+TEM10*TEM3)*TLFC 
TEM8= <  TLFA+TLFB*TEM5 ) *TEM5 
TEM9=TEM8*TLFP 
IF (TEM4.GT ♦ TEM9 ) TEM4-TEM9 
IF(TEM4«LT.-TEM9)TEM4=-TEM9 
TEM8=TEM8+TEM4 
IF  <  TEM8. GT ♦ TLFBR ) I SHU  < I )  =  1 
TPX-TPX+TEM8#TEM1 
TPY=TPY+T£M8*TEM2 
TPZ=TPZ+TEM8*TEM3 
G0T07500 
7450  I TEM=0 

C  AT  LEAST  ONE  SHROUD  LINE  IS  BROKEN 

D07455J=1 * NSHD 
7455  ITEM=ITEM+I3HD( J) 

IF ( ITEM«EG»NSHD)G0TC7550 
7500  CONTINUE 

C  TEAR  STRIP  BEHAVIOUR 

TEM1=SGRT<TPX*TPX+TPY#TPY+TPZ*TPZ> 

IF(TEM1«LT»  TLFT ) 60T07600 
IF(TLFL«GT. TLFTL ) G0T07530 
C  TEAR  STRIP  YIELDS 

TLFL=1 ,001*TLFL 
GOT07000 

7530  IFCTEhl .  L.T .  TLFF )  G0T07600 
C  PARACHUTE  HAS  BROKEN  FREE  FROM  BODY 

7550  TREL- 10000*0 
ISECON-l 
GOT05500 

C  TORSIONAL  TORQUE 

76 00  TEM2=TWA 

IF  ( TEM2«  GT  »3«0)TEM2-3.0 
IF  (  TEM2 . LT . -3 . 0 ) TEM2=-3  *  0 
TEM1=TEM1*TEM2*TWIST*SCLPAR*SCLPAR 
C  ADD  SHROUD  LINE  FORCES  TO  PARACHUTE  AND  BODY 

ADMPAR (1*7) = ADMPAR ( 1  *  7 )  +  1 PX 
ADMF'AR  (2*7)  -  ADMPAR  <  2 1 7 ) +TPY 
ADMPAR (3*7 )= ADMPAR (3*7) +  TPZ 
ADMPAR (4*7) -ADMPAR (4*7 ) +TEM1 
ADMPAR(5*7)=ADMPAR(5*7)-TPZ#TLX+TPX*TLZ 
ADMPAR (6*7) -ADMPAR  <6  *  7 ) +TPY#TLX-TPX#TLY 
TBX=-TPX*S1-TPY#S4-TPZ*S7 


TBY=*-TPX*S2-'TPY*S5-TF'Z*S8 

1  BZ»-TPX*S3-TPY*S6-TF'Z*S9 

ADMBOD  (1 »  7 )  bADMBOD  (1*7)  +TBX 

ADHBOD (2*7) -ADMBOD  (2*7) +TBY 

A  DM  BO  Ei  (  3 »  7  )  -ADMBOD  <3/7 )  +TBZ 

ADMBOD (4*7) =ADMBOD <  4 »  7 ) -TEM1*S1 

ADMBOD (5*7) =ADMBOD  <  5  *  7 > -TBZ*BODTL-TEMl *32 

ADMBOD ( 6 * 7 ) =ADMB0D<6*  7)+TBY*B0DTL-TEMl*S3 

CALL  ENERT  <  ADMBOD »  BDMBOD ,  RBOD ,  UBOD  *  DBOD ,  XBOD » YBOD  *  ZBOD  *  F'MBOD  * 
lF'YBOD  *  PZBOD  *  PI  XBOD  *  PI  YBOD  *  DT  ) 

CALL  ENERT  <  ADMF’AR  *  BDMPAR  *  RPAR ,  OPAR ,  DPAR  *  XPAR  *  YF‘AR  *  ZPAR  *  PMPAR  * 
10.0*0.0*  F'lXF'AR  *  F‘I  YF'AR » DT  ) 

CALL  DERI V  <  ADMBOD  *  OBOD  *  DBOD  *  XBOD  *  YBOD  *  ZBOD  *  DT ) 

CALL  DERI 0  ( ADMPAR  *  OPAR  *  DPAR  ,  XPAR  *  YF'AR  *  ZPAR  >  DT  ) 

C  PARACHUTE  DEPLOYMENT 

SCLPAR*SCLPAR+SDLPAR*<RF'AR<  1  >-ABS<RPAR<2)  > -ABS < RPAR < 3 )  )  >#DT 
TWA-TWA+  <  DBOD  <  4  ) -OPAR  <  4  )  )  *DT 
I F  <  SCLF'AR  tGT .  SCLMAX  )SCLPAR=SCLMAX 
IF ( 3CLPAR . LT . SCLMIN  ) SCLPAR=SCLMIN 
IF ( SCLOLD  .  EQ «  3CLMAX « AND  *  SCLPAR  .EQ.SCLMAX) 60T05700 
6005  TEM2=SCLPAR*SCLPAR 
D06100I=1 *  16 
PA  PAR  <  I )  =F'APBR  <  I )  #TEM2 
6100  F'WPAR  ( I )  -F'WF'BR  ( I )  *TEM2 
SCLOLD=SCLPAR 
F'ARRAD=SCLPAR*PBRRAD 
GO  T 05700 

C  BODY  ALONE 

5500  CALL  ENERT  < ADMBOD * BDMBOD  *  RBOD  *  DBOD *  DBOD  *  XBOD  *  YBOD . ZBOD  *  PMBOD  * 
1PYB0D  •  PZBOD . PI XBOD  *  PI YBOD  *  DT ) 

1 F'  (  IF IRST  4  EQ .  0  )  GOT 05505 

C  FIRST  PASS  THROUGH  PROGRAMME  TO  SET  UP  BDMBOD  ETC. 

IFl'RST-0 

G0T05705 

5505  CALL  DERI V < ADMBOD  » UBQD >  DBOD  *  XBOD  *  YBOD  *  ZBOD  *  DT ) 

C  INCREMENT  TIME 

5700  T  =  T  +  DT 

C  CHECK  PARACHUTE  STATUS 

IF  <  T . GT  . TREL . AND . ISECON . EQ , 1 ) G0T06000 
IF  <  T . GT . TREL1 ) SCLMAX=SCL2 
IF  <  T  .  GT  ♦  TREL.2  )  SCLMAX=  1 , 0 
I F  <  T  .  LT  *  PGI T ) G0T0200 
C  OUTPUT  A  RECORD 

C  THE  WRITE  STATEMENTS  AND  NEED  TO  CALL  EULER  SHOULD  BE 

C  OAR I  ED  AS  REQUIRED 

5705  PQIT-T+POIDT 

CALL  EULER ( TEM1 *  TEM2 *  TEM3 *  XBOD  *  YBOD  *  ZBOD ) 

IF < ISECON ♦ EQ . 1 >  GQT05900 

CALL  EULER  <  TEM4  *  TEM5  *  TEM6  *  XPAR  *  YF'AR  *  ZPAR ) 

WRITE  <  3  *  1 99  )  T  *  UBQD <1  )  *  OF'AR  <1 )  .  TEM2  *  TEM5  *  MBOD  (ID*  OBOD (12)*  UBOD  <  1 3  ) 
WRITE  <  7  *  199  ) T  *  VBQD<  1  )  * OF’AR<  1  )  * TEM2 *  TEM5*  UBOD (11)* UBQD ( 12  )  *  UBOD  ( 13  ) 
G0T05950 

5900  WRI TF.  (  3  *  1 99  )  T  *  OBOD  <  1  )  ,  0 . 0  *  TEM2  *  0 . 0  *  OBOD  ( 1 1  )  *  OBOD  <  12  )  *  OBOD  <  1 3  ) 

WR I TE  (  7  *  1 99  )  T  ,  OBOD  < 1  )  *  0 . 0  *  TEM2  *  0 , 0  *  OBOD  (ID*  OBOD  (12)*  OBOD  (13) 
FORMAT  <  F9  *  4  *  7F‘10  *  3) 


199 
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5950  IF <  T . LT  *  TMAX ) GOT 0200 
C  RUN  COMPLETE 

WRITE ( 3  *  199 ) -1 ♦ 0»  0 ♦ 0» 0 . 0 f 0 • 0 f  0 « 0»  0 ♦ 0»  0 *  0 » 0  « 0 

STOP 

END 

SUBROUTINE  XYZ(UfXfYfZ) 

C  GENERATES  THE  TRANSFORMATION  MATRIX  XfYfZ  FROM  QUATERNIONS 

DIMENSION  U<13) fX<3) f Y(3) fZ(3) 

X(1 )  =  1 ,0-2.0*<U<9)*V<9)+Vd0>*Vd0> ) 
X<2>=2.0*<Vd0>*V<7)+U<8)*0<9>  > 

X < 3 ) =2 . 0* ( V < 8 ) *U ( 1 0 > -0 < 7 ) *U < 9 ) ) 

Y<1 )=2.0*<U<8)*V<9)“V<7)*V<10> ) 

Y(2)  =  l  ♦0-2.0*<U(8)*V<8>+U<  10>*Ud0>  ) 
Y(3)=2.0*(M(7)*V(8)+y(9)*y(10)) 

Z<1  )=2.0*<V<7>*V<9)+V<8>*0d0)  ) 

Z<2)-2«0*<V<9>*VdO>-V<7>*V<8>  ) 

Z<3)-l«0-2«  0* (0<8)#V<8) +0(9 ) ( 9 ) ) 

RETURN 

END 

SUBROUTINE  SEA ( W f X f Y f Z f XPOS f YPOS f ZPOS ) 

C  MOTION  AND  GEOMETRY  OF  THE  SEA 

DIMENSION  W(6) f X(3) f Y<3) f Z<3) 

COMMON/BLOCK 1 /  T  f  SEAA f  SE AW  f SEAC  f  COMX  f  COMY  f  WDEP  f  SURMY  f  SURMZ  f 
1SQTEM1  f  ATSYSZf  SURECGf  IiURECGf  BETIDf  BET2Df  CGORf  XAXORf  CAUAf  CAVBf 
2SQVYUZf ATUZVYf BETlCfBET2Cf BOFUf BOFUfBOFWfPCAUf XAXfRADfDAXfSINAf 
3C0SA f I SW 1 f I SW2  f I TES 1 1 CROS  f I DEP  f I DRY  f  BODU  f  BODV  f  BODW  f  PROT 
ANG=SEAW*  <  XP0S*C0MX+ YPQS*COMY-SEAC*T ) 

S=SIN<ANG) 

C=COS  <  ANG ) 

ZEXP=SEAA*S 

I F  ( ZPOS ♦ GT ♦ ZEXP ) ZEXP=ZPOS 
E=£XP(-SEAW*ZEXP> 

A1=SEAA*SEAC*SEAW 

A2=A1*SEAC*SEAW 

C  SEA  MOTION  AT  XPOS f YPOS t ZPOS 

U=~A1*E*S 
Wd)=U*COMX 
W<2)=U*C0MY 
W(3)=-Al*E*C 
UD0T=A2*E*C 
W  <  4  )  -"UDOT&COMX 
W(5)-UD0T*C0MY 
W(6)=-A2*E*S 
WDEP-10 . 0+ZF'0S-SEAA*E*S 
C  SEA  SURFACE  GEOMETRY 

DZDX=SEAA*S£AW*C 
A1=SQRT(1.0+DZDX*DZDX) 

A2=DZDX/A1 
XN=A2*C0MX 
YN=A2*CQMY 
ZN=-1 • O/Al 

SURMY=XN*Y<1 )+YN*Y(2)+ZN*Y<3> 

SURMZ=XN*Z<1 >+YN*Z<2)+ZN*Z<3) 

TEM1=SURMY*SURMY+SURMZ*SURMZ 
IF(TEM1 .LT. 0.00001 )TEM1 =0,00001 
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SGTEM1=SQRT<TEM1) 

ATSYSZ=ATAN2 ( SURMY » SURMZ ) 

SURECG=ZN*<SEAA*S-ZPOS> 

DURECG=XN*X<1)+YN*X<2)4ZN*X<3> 

RETURN 

END 

SUBROUTINE  HYDRO <A»P> V*R) 

C  FULLY  IMMERSED  FORCES  AND  ADDED  MASSES 

DIMENSION  P(16)*A<6»7)»V<13)»R(6) 

U=ABS<R< 1 ) ) 

U2=R<1)*U 
UF'=U  <  4 )  *U 
UQ=U*V<5> 

UR=U*V<A> 

UY=U*R<2) 

UW=U*R<3) 

VWAB=SQRT  <  R ( 2 ) *R  <  2 ) 4R  <  3  >  *R  <  3  > ) 

U2=R<2)*UWAB 
U2=R<3)#VWAB 
GA2=U<5)*VWAB 
RA2=U<6)*VUAB 
C  FORCES 

A< 1 f 7 ) -P  < 1 ) *U24P <  1 4  >  #R  <  4  > 

A  ( 2  f  7 )  =P  <  3  >  *UR4P  <  2 >  *UV4P  <1 0  )  *U24F*  <  1 1 >  *RA2+P  ( 1 4  )  *R  ( 5  ) 

A(  3 » 7 )  =~P  < 3 )  #UQ4F‘  ( 2 ) #UW4F*<  10) #W2-F‘<  1 1 ) #GA24F* < 14 ) #R  <  6 ) 

A<4»7)=P< 1A)*UP 

A( 5  f  7 )  =  -P < 1 5 ) *R  <  6 ) +P ( 4 ) *UW4P < 5 ) *UG4P  <12 ) #W2+P (13) *GA2 
A<  6 , 7 )  «F  <  1 5 )  *R  <  5 )  -P  <  4  )  *UV4P  <  5  >  *UR-F*  <12 )  *V24P  <  13  )  *RA2 
C  ADDED  MASSES 

A<1»1)=P<A) 

A<2>2 )-P<7 ) 

A<  2  r  6  ) -F  <8 ) 

A< 3  > 3  )  =F'  <  7 ) 

A<  3 » S ) =-P  <  8 ) 

A<5»5)=F'<9) 

A<  A » A  )  -F' < 9 ) 

RETURN 

END 

SUBROUTINE  DEPTH 

C  CALCULATES  THE  INTERSECTION  OF  A  SEGMENT  WITH  THE  SEA  SURFACE 

C  BETID  TO  BET 2D  IS  THE  INTERSECTION  SECTOR 

COMMON/BLOCK 1/  T »  SEAA*  SEAW»  SEAC»COMX  » COMY  t WDEF' t SURMY  r SURMZ » 
1S0TEM1 »  ATSYSZ  » SURECG »  DURECG » BET  ID  >  BET2D » CGOR  >  XAXOR »  CAUA » CAVB  t 
2SQUYVZ  f  ATOZOY  t  BET  1C » BET2C»  BOFU  r  BQFV  »BQFW » F’CAV » XAX » RAD»  DAX » SINA » 
3C0SA  » ISW1 * ISW2  > ITES»  ICROSr  IDEF'r  I  DRY  t  BQDU  r  BODY » BQDW»  F'ROT 
IDEF'-O 
I  DRY-0 

SURE«SURECG-XAX*DURECG 

SURE-SURE/SGTEM1 

IF  <  SURE ♦ LE  *  RAD ) GOTO 120 

XDEP-1 

RETURN 

120  IF ( SURE . LE . -RAD) GOT 0500 
SURE=SURE/RAD 
TEM1 =SQRT ( 1 . 0-SURE*SURE ) 
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T£M2=ATAN2 ( SURE  t TEM1 ) 

B£T1D=TEM2~ATSYSZ 
B£T2D=3, 1416-TEM2-ATSYSZ 
TEMl=<BETlD+BET2D>/2.0 
TEM2=C0S<TEM1)*SURMY+SIN<TEM1)*SURMZ 
IF ( TEM2.LT. 0.0)RETURN 
TEM1=BET1D 
BETiIi*BET2B 
BET2D=TEM1 
RETURN 
500  IDRY*1 
RETURN 
END 

SUBROUTINE  FORCES ( A  *  R » HT , CL IFT , CBRAG ) 

C  FORCES  AND  ADDED  MASSES  FOR  EACH  WET  SECTOR 

DIMENSION  A(6»?),R<*> 

COMMON/BLOCK  1  /  T » SEAA »  SEAW » SEAC  t  COMX  >  COMY » WDEF' » SURM Y » SURMZ  f 
1SGTEM1  t  ATSYSZ » SURECG*  DURECG*  BET1D»  BET2D»  CGOR  t  XAX0R»  CAVA » CAVB» 
2SGVYVZ  r  ATVZVY  t  BET  1C  ?  BET2C » BQFU»  BQFV » BOFW  t  PCAV»  XAX » RAD » DAX » SINA  t 
3C0SA 1 1SW1 1  ISW2»  ITESf  I CROS  » IDEF* »  IDRY >  BODU ,  BODV » BODW »  PROT 
C  SORT  OUT  ACTUAL  WETTED  SECTOR  FROM  DEPTH  AND  CAVITATION  DATA 

I SEG-0 

IF  < IDEP*EG« 1 )G0T0300 
IF(ITES.EG*0) G0T0205 
BET1C=BET1D 
BET2C«BET2D 
G0T0300 

205  TEM2=SCALE  C  BET2D-BET1D ) 

TEM3«SCALE  <  BET 1C-BET 1 D ) 

TEM4=SCALE  <  BET2C-BET1D) 

220  IF(T£M4.LT,TEM3)G0T0240 

IF  ( TEM3 « GT  *  TEM2) RETURN 
BET1C=TEM3+BET1D 
IF ( TEM4 • GT  »  TEM2) G0T023O 
225  BET2C=TEM4+BET1D 
GOT0300 

230  BET2C=TEM2+BET1D 
GOT0300 

240  BET1C-BET1D 

IFCTEM4.GT «  TEM2 ) GQT0230 
IF ( TEM3 . GT . TEM2 )G0T0225 
BET2C=TEM4+BET1D 
ISEG=1 

BET3C=TEM3+BET1D 

BET4C=TEM2+BET1D 

C  CALCULATE  FORCES  AND  ADDED  MASSES 

300  BET l C=SC ALE (BET 1C ) 

BET2C=SCALE ( BET2C ) 

IF( BET2C «LT* BE TIC) BET2C=BET2C+6« 283185 
SINB1=SIN<BET1C) 

SINB2=SIN (BET2C ) 

C0SB1=C0S(BET1C) 

C0SB2=C0S ( BET2C ) 

IF(ISW1«EG«1*  AND  * ISW2 « EQ. 0) G0T0500 
C  AI  TO  FI  INTEGRALS  USED  FOR  FORCES  AND  MASSES 


AI=BET2C-BET1C 
BI-SI NB2-SINB1 
CI=C0SB1 -C0SB2 

TEM3= (SIN ( BET2C+BET2C ) -SIN <  BET 1C+BET1C ) >  /4  4  0 

TEM4=AI/2.0 

TEM7=SIN<TEM4) 

IF(AI.GT.3, 142)TEM7=1 ,0 
DI=TEM4+TEM3 

£I=<SINB2*SINB2-SINB1#SINB1 )/2,0 

FI=TEM4-TEM3 

TEM1=RAD*DAX 

TEM2=1030,0*T£M1 

TEM3=ABS ( COSA ) 

TEM4=SINA/TEM3 
TEM5=C0SA/TEM3 
TEM6=RA0*SINA-XAX*C0SA 
TEM7-TEM2*HT*TEM7 
TEM8=TEM7*TEM4 
TEM9=TEM8*SINA 
TEM10=TEM8*COSA 
TEM11=TEM8*TEM6 
TEril2-TEh7*TEM5 
TEM13=TEM12*C0SA 
TEM14=TEM12#TEM6 
TEM1S=TEM7*TEM6*TEM6/TEM3 
ADDED  MASSES 
ACL  1 1  >=A<  1 1 1  )+TEM9#AI 
A(  1 f2)=A( 1 >2)+TEM10*BI 
A<1»3)=A<1»3) 4TEM10*CI 
A<1»5>=A<1>S)4TEM11*CI 
AU  >6)=A(  1  »6)-TEMll*BI 
A  <  2 » 2 ) =A  <  2  >  2 ) 4T£M13*DI 
A ( 2  » 3 ) =A ( 2 »  3 ) 4TEM13*E  I 
A<2?5)-As2*5)+T£M14#EI 
A ( 2 » 6 ) =A  <  2  >  6 ) -TEM14KDI 
A<3>3)=A<3,3)4T£M13*FI 
A<3»5)=A(3»5)4TEM14*FI 
A<S,5)=A<5>5>4T£M15*FI 
A(Sf6)~A<5f 6>-TEM15*EI 
A(6»6)-A<A»6)+TEM15*DI 
TEM8- <  BET 1C4BET2C ) *0  4  5 
TEM7---C0S  <  TEM8 ) 

TEM8=SIN<TEM8) 

TEM7-BOFU40 . 25* <  BOFV* ( C0SB1 4C0SB24TEM74TEM7 ) 4BQFW* 
1  ( SINB14SINB24TEM84TEM8) ) 

TEM8~CDRAG*TEM2*ABS(TEM7> 

TEM9-TEM8#TEM5 

TEM10=B0FU*AI4B0FV*BI4B0FW#CI 
TEM11~B0RJ*BI4B0FV*DI4B0FW*EI 
TEM12=B0FU*CI4B0FV*EI4B0FW*FI 
I F ( I SW2 « EQ  4  0 ) G0TG350 
NON  LINEAR  FORCE 

DAFX~TEM8*TEM4* ( TEM10+T£M10-TEM7*AI ) 

DAFY®TEM9* (TEM114TEM11 -TEM7*BI > 
DAFZ*TEh9#(TEhl2+TEM12-TEM7*CI ) 


IF( ISW1 .EG* 1 )GQTQ400 

C  BUOYANCY  FORCE 

TEM11=0.5#TEM2*RAD*<AJ-SIN<AI ) ) 

DAFX=DAFX-R<4)*TEMU 
OAF Y=DAFY-R ( 5 ) #TEM1 1 
DAFZ=DAFZ-R<6)*TEM1 1 
GQT0400 

C  LINEAR  FORCE 

350  TEM7=CLIFT*TEM2*SGRT  <  BODU*BODU+BODU*BOnU+BODW#BODW ) 

TEMS=TEM7*TEM5 
DAFX=TEM7#TEM4*TEM10 
DAFY=TEM8#TEM11 
DAFZ=TEM8*T£M12 

400  A( 1 »7)=A( 1 f 7)-DAFX 
A(2»7)=A<2»7) -DAFY 
A<3f7)=A<3»7) -DAFZ 
TEM14=TEM6/C0SA 
A(5f7)=A(5»7)-DAFZ*TEM14 
A(6»7)=A<6»7)+DAFY*TEM14 
G0T0900 

C  CRUCIFORM  TAIL 

500  TEM6=1030.0*DAX 

TEM10=TEM6*CLIFT*SGRT<BQOU*BOnU+BQBU*BOnU+BOnW#BQDW 
TEM 1 =RAD#RAD#0 ♦ 5 
T£M2=T-0 ♦ 5~FL0AT  <IFIX<T) ) 

TEM3=TEM1*0,01 

TEM4=0. 3/RAD 

TEM5=0. 4/RAD 

Y1--RAD 

Y2=RAD 

Z1=-RAD 

Z2=RAD 

IF(SINB1#SINB2.GE.0.0) G0T0520 

TEM5=RAD*<SINB1*C0SB2-CQSB1#SINB2)/<SINB1-SINB2> 

IF(SINB1.GT«0.0)G0T0510 

Y1=TEM5 

G0T0540 

510  Y2=TEM5 
G0T0540 

520  IF<  < SINB1 +SINB2 ) . GT . 0 . 0 )GOT0530 

IF< (BET2C-BET1C) , GT . 3 . 1 4 1 59 > G0T0540 
G0T0700 

530  IF  <  C0SB1 ♦ GT ♦ CQSB2 ) GOTO 700 

540  AI= Y2-Y 1 

BI=AI*(Y2+Y1 )/2,0 
CI=<Y2*Y2*Y2-Y1*Y1*Y1  )/3.0 
T£M7=TEM6#HT*AI/ ( RAD+RAD ) 

TEM8=TEM7#XAX 

C  ADDED  MASSES  AND  FORCES  'HORIZONTAL  FIN' 

A(3»3)=A(3»3)+TEM7*AI 
A<3>4)=A<3>4)+TEM7*BI 
A<3»5)=A(3>5)-TEM8*AI 
A(4»4)=A<4>4)+TEM7#CI 
A(4»5)=A(4»5)-TEM8*BI 
A<5>5>=A<5»5)+TEM8*XAX*AI 


J4. 

DAFZ=TEM10*(B0DW*AI+F‘R0T*BI ) 

A(3»7)SSA(3»7)  -EiAFZ 

A(4»7)SSA<4»7)-TEM10*( BQDW#BI+F'ROT#CI ) 

IF  <  ABS( BI > .LT*TEM3)A<4>7)=A<4»7)+ 

1<SIGN<TEM4>  TEM2 ) -PROT  ) #ABS  < DAFZ  )  HcTEMl 
A<5t7)=A<5>7)+XAX*DAFZ 
700  IF ( C0SB1 *C0SB2 . GE . 0 , 0  > G0TG720 

T£M5=RAD* ( S INB2*C0SB 1 -C0SB2*SINBI ) / ( COSB 1 -C0SB2 ) 

IF<C0SB2*GT ♦0*0)GGT0710 
Z 1 =TEM5 
GCT0740 
710  Z2=TEK5 
G0T0740 

720  IF ( <  C0SB1 +C0SB2 ) ♦ GT ♦  0 «  0 ) GOTO 730 

IF< ( BET2C-BET1C ) «GT *3.14159) GO TO 740 
G0T0900 

730  IF  <SINB2«GT .SINB1 )GGT0900 
740  AI=Z2-Z1 

BI=AI*<Z2+Zl)/2.0 
CI=(Z2*Z2*Z2-Z1*Z1*Z1 )/3.0 
TEM7~T£M6*HT*AI/<  RAD+RAD) 

TEM8=TEM7*XAX 

C  ADDED  MASSES  AND  FORCFJS  'VERTICAL  FIN" 

A<2»2)=A<2>2)+TEM7*AI 
A<2»4)=A(2,4)-TEM7*BI 
A<2»6)=A(2>6 ) +T£M8#AI 
A(4f4)=A(4,4)+TEM7*CI 
A<4>6)=A<4,6)-TEM8*BI 
A < 6  *  6 ) -A < 6 , o ) +TEM8*XAX*A I 
DAFY-TEM10*<BGDV*AI-PR0T*BI ) 

A  <  2  f  7  )  =A  <  2  •  7 )  -DAF'Y 

A(4>?)=A<4r7)+TEM10*<BGDV*BI-F'R0T*CI) 

ZF<AE<S<BI ) •  LT .TEM3)A(4*7)=A(4»7>+ 

1 <  SI GN  <  TEM5 »  TEM2 ) - PROT ) #ABS <  DAFY  >  *TEM1 
A(6»7)=A(6>7) -XAX#DAFY 
900  IF< ISEG»EQ«0 i RETURN 
BET1C-BET3C 
BET2C=BET4C 
ISEG-0 
G0T0300 
END 

SUBROUT  I NE  ENERT  (  A  t  B »  R  >  0 > D  >  X  r  Y  >  Z » PM  f  PY ,  F'Z  » P IX » F* I Y » DT  ) 

C  FORCES  DUE  TO  GRAVITATIONAL  AND  INERTIAL  SOURCES  BUT 

C  EXCLUDING  TERMS  CONTAINING  DU/DT»DV/DT» . . DR/DT 

DIMENSION  V(13)»X(3)fY(3)»Z(3)»D<6)»S(6)fR(6)»A(6>7)»B(6r6) 
F2=V<4)*V<4) 

F'G=V(4)*V<5> 

PR=-'V<4)*V<6) 

G2=V<5)*V<D) 

R2-V ( 6 ) *V  <  6 ) 

QR-V<5)*V<6) 

UQ-V< 1 >*V<5> 

UR-V(6)*V(1 ) 

VP“V(4)*Vv2> 

VR=V<2)*V<6> 


35. 


WP<*V<3>*V(4) 

WG=V<3>*V<5) 

YM=PM*PY 
ZM»PM*PZ 
X3=X<3>*9,81 
Y3=Y < 3) *9 <  81 
Z3=Z<3>*9,81 
C  MOTION  FORCES 

A  <  1 » 7 ) =A  < 1 >  7 ) -PM* ( WG-VR+PY*PG+PZ*PR~X3 ) 

A<2f7)»A<2»7>  -PM*  <  UR-Wr'-PY*  ( R2+P2 )  +PZ*GR- Y3 ) 
A(3f7)=A<3»7)  -PM*  ( VP-UG+PY*GR-PZ*(F'2+Q2 )  -Z3 ) 
A<4»7)=A<4*7)-ZM*( WP-UR+Y3) -YM*< VP-UQ-Z3 ) 
A(5»7)=A(5»7)-(PIX-PIY)*PR-ZM*(WQ-VR-X3) 

A ( 6 , 7 ) =A ( 6 ,  7 ) - <  P I Y-P I X  > *PG-YM* <  UR-WQ+X3 ) 

C  ADDED  MASS  EQUALITIES 

A<2»1>®A(1»2) 

A(3r 1 )=A( 1 »3) 

A(3»2)t=A(2»3) 

A(3»6)=-A<2»5) 

A ( 4  » 2) =A ( 2 »  4 ) 

A(4r3)=A(3*4) 

A<5*1)=A(1,5) 

A  <  5 » 2 ) =A<  2 »  5 ) 

A(5»3)=A<3*5) 

A<5»4)=A(4»5> 

A(6»l)=A(li6) 

A ( 6 » 2 ) “A <  2 » 6 ) 

A  ( 6  >  3 ) =A( 3 »  6 ) 

A (6 » 4 ) =A (4*6) 

A<6,5)=A(5 ,6) 

S < 1 ) =R<  4 ) +X3 
S < 2 ) =R<  5 ) + Y3 
S < 3 ) =R<  6 ) +Z3 
D03I=4 1 6 
S<I)=0.0 

3  R( I )=U< I ) 

TEM1=A  < 1 » 1 )+A<2>2)+A<3>3) 

TEM2=B< 1 »1 )+B<2f 2)+B<3»3> 

IF  <  TEM1 • GE . TEM2 ) GQT04 
C  SHEDDING  MASS 

D0100I=1 f 6 
D0100J=1 » 6 

100  A<I» J)=B( I » J)+0.2*< A< I f J)-B< I f J) > 

4  DOS 1  =  1 »  6 

3  D< I )=0*  0 

00101  =  1 » 6 

C  ADDED  MASS  FORCES 

A  < 1 » 7 )=A( 1 » 7  >-R< I ) * ( A<3 » I >*R( 5 ) -A  <2» I >*R  <  6 ) )+A( 1 » I ) *S  1 1 ) 
A<2»7)-A(2»7)-R<I)*<A<1»I)*R(6)-A(3»I)*R(4))+A<2»I)*S(I> 
A ( 3 » 7 ) =A ( 3 f  7 ) -R  < I ) * ( A  <  2  *  I ) *R ( 4 ) - A  < 1 f I ) *R ( 5 ) ) +A  <  3 » I ) *S ( I ) 
A ( 4 » 7 ) =A ( 4 » 7 ) -R ( I ) # ( A ( 3 » I ) *R ( 2 ) -A  <  2  f I ) *R  <  3 ) +A  <  6 1 1 ) *R  <  5  > 
1-A  ( 5  >  I )  *R  ( 6  )  )  +A  <  4 » I  )  *S  (  I  > 

A ( 5 » 7 ) =A ( 5 »  7 ) -R  < I ) * ( A  < 1 » I >  *R ( 3 ) -A  <  3 » I ) *R ( 1 ) +A ( 4  f I ) *R  <  6 ) 
l-A(6fI)*R(4) )+A(5fI )*S(I) 

A<6f7)=A<6»7)-R<I)*<A<2»I)*R(l>-A<l»I)*R<2)+A(5»I>*R<4) 
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l-A<4rI)*R<5> >+A<6,I>*S<I> 

D010J--1  *& 

IF  <  TEM1 « LT . TEM2 ) GOTO 10 
C  RATE  OF  CHANGE  OF  MASS 

D<J)=D<J>-R<I>*<A<J»I>-B<JrI>> 

C  OLD  ADDED  MASSES 

10  B( J» I )=A( J» I ) 

C  TOTAL  MASS  MATRIX 

Aa»i>=A<l»l)+PM 
A(li5)=A(l»5)+ZM 
A<1.6)=A<1>6)-YM 
A(2*2)=A(2»2) +PM 
A  (  2  >  4  ) -A <2*4 ) -ZM 
A  ( 3  *  3 ) -A (3*3) +PM 
A(3,4)=A<3*4)+YM 
AC4»2)=A(4»2)-ZM 
A  ( 4  *  3 )=A<  4 *  3 ) +YM 
A  (  4  »  4  )=3A<  4  *  4 ) +PI X 
A  ( 5  *  5 ) =A<  5  *  5 ) +PI Y 
A(5, 1 )=A<5* 1 >+ZM 
AC6*1)=A<6, 1)-YM 
A(6»6)=A(6r6)+PIY 

C  ADJUST  TIME  INCREMENT  TO  SUIT  NATURE  OF  MOTION 

TEM1=ABS(D<1)/A(1»1))+ABS<D<2)/A(2»2)) +ADS (D(3)/A(3»3) ) 
1+50.0*<ABS<DC4)/A<4*4) >+ABS<D<5>/A<5*5> ) +ABS ( D ( 6 ) /A ( 6  *  6 ) ) ) 

C  SET  DT=0  FOR  IMPULSIVE  VELOCITY  CHANGE 

IF<TEM1.GT.1.0)DT=0.0 

TEM1~ABS(A< 1 »7)/A(l» 1 ) ) +ABS < A< 2* 7 > /A ( 2 > 2 > >+ABS< A < 3 » 7 ) /A ( 3 » 3 ) ) 
1+50.0*(ABS(A(4» 7)/A<4»4) )+ABS<A<5*7)/A<5*5> ) +ABS ( A<6 * 7 ) /A < 6  *  6) > ) 
TEM1=0 «  05/SQRT ( TEM1 ) 

C  REDUCE  DT  IF  ACCELERATION  IS  LARGE 

IF(DT.GT.TEM1)DT=TEM1 
RETURN 
END 

SUBROUTINE  DER IV(A*V»D>X*Y»Z*DT) 

C  SOLVE  AND  INTEGRATE  THE  EQUATIONS  OF  MOTION 

DIMENSION  A<6f7)»V(13)»D(6)»X(3)»Y(3)»Z(3)tDV(6) 

D03I-1 ? 6 

C  EVALUATE  HALF  TOTAL  EXTERNAL  FORCE  IMPULSE 

3  A ( 1 1 7 ) =0 ♦ S*  <  A  < 1 1 7 ) *DT+D < I) ) 

C  SOLVE  EQUATIONS  OF  MOTION  BY  FORCING  THE  LOWER  L.H.  OF 

C  THE  MASS  MATRIX  TO  ZERO  AND  THE  DIAGONAL  TO  UNITY 

0020J=1  ,5 
TEM=A  <  J  >  J ) 

D010N= J * 7 

10  A( JfK)=A< JfK)/TEM 

I=J+1 
D02QL=I » 6 
D020K=I » 7 

20  A<L,K)=A<L*K)-A<Lr J)*A< J*N> 

DV ( 6 ) =A  <  6  *  7 ) /A (6*6) 

DV  <  5 ) =A  <  5  *  7 ) -A  <  5*6) *DV ( 6 ) 

DV(4)=A<4* 7)-A<4*6)*DV<6)-A<4*5)*DV(5> 

DV  <  3 ) “A  <  3  *  7  > -A ( 3 * 6 ) *DV <  6 ) ~A  <  3  *  5 ) *DV C  5  )  -  A  <  3  *  4  >  #DV <  4 ) 

DV ( 2 ) = A  <  2  *  7 ) -A ( 2 * 6 ) *DV ( 6 ) -A ( 2 * 5 ) *DV < 5 ) -A  <  2  *  4 ) *DV ( 4 ) -A  <  2 *  3 ) *DV <  3 ) 


DV  < 1 ) =A < 1 , 7 ) -A  <1 » 6 ) *DV  <  6 > -A  < i 1 5 ) *DV  <  5 ) -A  < 1  *  4 ) *BV (4 )-A< 1 »3) *DV  <  3  > - jfl 
1A< 1 »2)*DV<2)  1 

11058501  =  1 » 6 

C  INCREMENT  VELOCITIES  BY  HALF  TOTAL  INCREMENT 

5850  V  < I ) =V  < I >  +DV ( I > 

DT2=nT#0.5 

C  INCREMENT  QUATERNIONS 

V < 7 ) =V ( 7 ) - < V < 8 ) *V ( 4 ) +V < 9 ) *V < 5 ) +V < 10 ) *V < 6 ) ) *DT2 

V<8)=V<8)+<V<7)*V<4)-V< 10>*V<5>+V<9>#V<6> )#DT2  * 

V < 9 ) =V ( 9 ) + < V (10 ) *V < 4 ) +V < 7 ) *V < 5 ) -V < 8 ) *V ( 6 ) ) *DT2 

V  si 0 ) =V ( 1 0 ) - ( V < 9 ) * V ( 4 ) - V < 8 ) *V ( 5 ) - V < 7 ) #V ( 6 ) ) *DT2 

V  ( 1 1 ) =V < 1 1 )  +  < X ( 1 ) *V < 1 ) + Y ( 1 ) *V < 2 )  +  Z < 1 ) *V < 3 ) ) *DT 
V<12)=V<12)+<X<2)#V<1)+Y<2)*V<2)+Z(2)*V<3>  >*DT 

V(13)=V<13)+<X(3)#V<l)+Y(3)#V(2)+Z<3)#V(3))!fcDT  J 

11058601  =  1  >6  i. 

C  INCREMENT  VELOCITIES  BY  REMAINING  HALF  INCREMENT 

5860  V ( I >  =V  < I ) +DV ( I ) 

TEM=(3.0-V<7)*V<7)~V<8>*V<8>-V<9)#V<9)-V<10>*V<10>  >/2.0 
D05870I  =  7  > 10 

C  NORMALISE  QUATERNIONS  i\ 

5870  V<I)=V<I)#TEM 
RETURN 
END 

FUNCTION  SCALE ( X ) 

C  MANE  ANGLE  X  IN  RANGE  0  TO  2*F‘I 

:L0  IF<X.GE,0.0)G0T020 

X=X+6* 283185 

G0T010  | 

20  IF  <  X . LT . 6 ♦ 283185 ) 60T030 

X=X-6. 283185 
G0T020 

30  SCALE=X 

RETURN 
ENr' 

SUBROUTINE  EULER ( S » T > F » X » Y » Z ) 

C  EXTRACT  ROLL  <  F  )  >  F'ITCH<  T  >  > AZIMUTH  (  S )  IN  DEGREES  FROM  THE 

C  TRANSFORMATION  MATRIX 

DIMENSION  X<3) »Y<3) » Z(3)  ■ 

T1=X(3)*X<3) 

T2=0.0 

IF<T1 .LT. 1 .0)T2=SQRT< 1 .O-Tl ) 

T=~57 ♦ 3#ATAN2 ( X  <  3 ) r  T2 ) 

F=57*3*ATAN2(Y<3> »Z(3) ) 

S=57.3*ATAN2<X<2) ,X<1 > >  , 

RETURN 

END 
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FIG.  2.  A  TYPICAL  LIGHT  WEIGHT  TORPEDO 
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